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Abstract 

An  overview  of  current  unstructured  mesh  generation  and  adaptivity  techniques  is  given. 
Basic  building  blocks  taken  from  the  field  of  computational  geometry  are  first  described. 
Various  practical  mesh  generation  techniques  based  on  these  algorithms  are  then  constructed 
and  illustrated  with  examples.  Issues  of  adaptive  meshing  and  stretched  mesh  generation 
for  anisotropic  problems  are  treated  in  subsequent  sections.  The  presentation  is  organized 
in  an  educational  manner,  for  readers  familiar  with  computational  fluid  dynamics,  wishing 
to  learn  more  about  current  unstructured  mesh  techniques. 
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1  Introduction 


The  generation  of  unstructured  meshes  for  computational  fluid  dynamics  problems  has  evolved 
rapidly  over  the  last  ten  years.  In  fact,  it  has  been  stated  recently  that  unstructured  mesh  generation 
has  reached  such  a  level  of  maturity  that  it  can  be  considered  a  solved  problem.  Certainly,  the 
generation  of  unstructured  meshes  about  arbitrarily  complex  three-dimensional  configurations  can 
be  routinely  performed  on  present-day  workstations,  with  currently  available  techniques.  In  fact, 
a  reader  with  good  programming  skills  and  persistence  should  be  able  to  create  his/her  own  mesh 
generation  code  from  the  information  contained  in  these  notes.  However,  it  would  be  premature  to 
state  that  further  improvements  are  not  needed  in  this  field.  Of  the  various  techniques  described  in 
these  notes,  each  has  its  particular  strengths  and  weaknesses.  Further  improvements  in  efficiency, 
and  especially  robustness  are  surely  needed.  A  good  definition  of  an  optimal  mesh,  as  this  relates 
to  numerical  properties  of  the  solution  scheme  is  still  lacking,  and  with  it,  appropriate  techniques 
for  constructing  such  optimal  meshes,  particularly  for  anisotropic  or  stretched-mesh  generation. 

Many  of  the  successful  algorithms  in  unstructured  mesh  generation  have  found  their  roots  in 
the  field  of  computational  geometry.  Computational  geometry  is  the  theoretical  science  concerned 
with  defining  or  postulating  the  existence  of  specific  geometrical  constructs  [i.e.,  particular  tri¬ 
angulations  in  our  case),  devising  algorithms  for  generating  these  constructs,  and  analyzing  the 
complexity  of  these  algorithms  (usually  asymptotic  worst  case  complexities).  For  example,  the 
Delaunay  triangulation  represents  a  fundamental  computational  geometrical  construct  for  which 
many  construction  algorithms  have  been  devised  and  analyzed  [1.  2].  Unfortunately,  much  of  the 
computational  geometry  work  has  been  confined  to  two-dimensional  triangulations,  which  consti¬ 
tute  planar  graphs,  and  are  thus  easier  to  analyze.  Also,  the  concept  of  an  optimal  triangulation 
or  algorithm  from  a  computational  geometry  point  of  view  does  not  always  coincide  with  the  view 
from  a  computational-fluid-dynamics  point  of  view,  and  thus  many  of  the  computational  geom¬ 
etry  results  have  found  little  use  in  the  area  of  unstructured  mesh  generation.  Great  progress 
in  unstructured  mesh  generation  has  been  made  by  devising  heuristic  algorithms  combined  with 
empirical  experience.  In  some  sense,  the  engineering  field  of  mesh  generation  has  outstriped  the 
more  theoretical  field  of  computational  geometry,  particularly  for  three-dimensional  constructions. 
However,  while  heuristic  algorithms  may  work  well  for  most  cases  encountered  in  practice,  the 
lack  of  any  sound  theoretical  justification  of  these  algorithms  leaves  the  door  open  for  possible 
situations  which  may  result  in  failure  or  greatly  increased  complexity  of  the  algorithm.  Hence,  the 
construction  of  efficient  and  particularly  robust  mesh  generation  techniques  can  only  be  achieved 
through  algorithms  with  solid  theoretical  justification. 

The  purpose  of  these  notes  is  to  familiarize  the  reader  with  current  techniques  for  unstruc¬ 
tured  mesh  generation,  while  exposing  their  individual  strengths  and  weaknesses.  A  preliminary 
discussion  of  selected  computational  geometry  algorithms  which  are  most  relevant  to  mesh  gen¬ 
eration  techniques  is  first  given.  Once  these  essentials  have  been  outlined,  a  critical  description 
of  various  mesh  generation  approaches  in  both  two-  and  three-dimensions  is  given.  The  issue  of 
mesh  adaptation  is  discussed  in  a  following  section.  In  the  final  section,  modifications  to  the  de¬ 
scribed  algorithms  as  well  as  alternate  strategies  for  generating  anisotropic  or  stretched  meshes  are 
described. 

The  data-structures  required  for  efficient  implementation  of  many  of  these  algorithms  are  not 
discussed  in  detail.  Their  description  and  analysis  can  be  found  in  appropriate  computer  science 
text  books.  The  main  data-structures  employed  in  the  algorithms  of  this  chapter  are  the  heap-list, 
and  the  region  quadtree.  The  heap-list  represents  a  particular  implementation  of  a  priority  queue, 
which  enables  the  ordering  of  elements  based  on  a  key  [.3].  The  region  quadtree,  or  octree  in  three 
dimensions,  is  a  spatial  decomposition  data-structure,  which  enables  efficient  implementation  of 
spatial  search  operations  such  as  proximity  searches.  Quadtrees  are  discussed  in  detail  in  [4].  Both 
the  quadtree  and  the  heap-list  are  dynamic  data-structures,  i.e.,  they  enable  simple  operations 
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such  as  insertion  and  deletion  of  elements,  and  are  thus  particularly  well  suited  for  mesh  generation 
and  adaptation  purposes. 


2  Computational  Geometry  Constructs  and  Algorithms 


2.1  The  Delaunay  Triangulation 


Given  a  set  of  points  {P}  in  a  plane,  there  exists  many  ways  to  join  the  points  together  to  form  a 
set  of  non-overlapping  triangles  which  completely  covers  the  domain.  The  Delaunay  triangulation 
represents  a  particular  construction  of  this  type  which  has  various  well  defined  properties.  For 
example,  the  Delaunay  triangulation  is  the  dual  of  the  Voronoi  Tessalation.  A  Voronoi  Tessalation 
is  the  graph  obtained  by  drawing  the  median  line-segments  which  separate  the  plane  into  regions 
which  are  closer  to  a  given  point  of  {P}  than  to  any  other  point  in  the  set  {P},  as  shown  in  Figure 
1. 


Figure  1:  Voronoi  tessalation  and  corresponding  Delaunay  triangulation  of  a  set  of  points 
in  the  plane,  illustrating  the  empty  circumcircle  property. 

If  we  draw  a  line  segment  between  any  two  points  which  are  neighbors  in  this  Voronoi  diagram,  we 
obtain  the  Delaunay  triangulation  of  these  points. 

Another  property  of  the  Delaunay  triangulation  is  known  as  the  empty  circumcircle  property. 
This  states  that  no  point  of  the  forming  set  {P}  can  be  contained  inside  the  circumcircle  of  any 
triangle,  as  shown  in  Figure  1. 

It  is  also  well  known  that  the  Delaunay  triangulation  corresponds  to  a  max-min  triangulation, 
i.e..  out  of  all  the  possible  triangulations  of  a  given  set  of  points,  it  is  the  triangulation  which 
incurs  the  largest  minimum  angle  for  all  triangular  elements.  Thus,  a  Delaunay  triangulation  may 
be  expected  to  result  in  “well-shaped”  elements,  without  very  small  angles. 

While  there  are  many  other  properties  associated  with  the  Delaunay  triangulation  [1,  2],  the 
above  are  the  most  useful  for  mesh  generation  purposes. 

The  empty  circumcircle  property,  in  particular,  is  attractive,  since  it  extends  readily  to  three 
dimensions,  by  considering  the  circumsphere  associated  with  each  tetrahedron.  The  empty  circum- 
circle/sphere  property  forms  the  basis  for  several  Delaunay  triangulation  algorithms  in  both  two 
and  three  dimensions,  which  are  described  below. 

2.1.1  Bowyer- Watson  Algorithm 

This  is  an  incremental  algorithm  which  assumes  new  points  are  added  sequentially  to  an  existing 
Delaunay  triangulation  [5,  6].  When  a  new  point  is  inserted  into  the  triangulation,  the  first  task 
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is  to  identify  all  existing  triangles  whose  circumcirde  is  intersected  by  the  new  point,  as  shown 
in  Figure  2.  This  can  be  achieved  by  first  locating  the  triangle  which  contains  the  new  point. 
The  circumcirde  of  this  triangle  must  certainly  be  intersected  by  the  new  point.  The  remaining 
intersected  triangles  can  be  identified  by  first  searching  the  neighbors  of  the  containing  triangle, 
and  then  the  neighbors  of  these  neighbors,  in  a  breadth-first  type  search,  abandoning  the  search 
along  any  neighbor  path  where  the  test  result  is  negative.  Properties  of  the  Delaunay  triangulation 
guarantee  that  such  a  neighbor  search  is  sufficient  to  locate  all  intersected  triangles.  The  union  of 
these  intersected  triangles  is  then  removed  from  the  current  triangulation,  thus  creating  a  convex 
cavity  which  contains  the  new  point,  as  shown  in  Figure  2.  A  new  triangulation  is  then  constructed 
by  joining  the  new  point  to  all  vertices  on  the  boundary  of  the  polygonal  cavity.  The  algorithm 
extends  naturally  to  three  dimensions  by  considering  the  circumspheres  of  tetrahedra  and  results 
in  the  retetrahedralization  of  a  convex  polyhedral  cavity.  Proofs  of  the  convexity  of  the  resulting 
cavity,  and  the  validity  of  the  retriangulation  of  this  cavitv  have  been  formulated  in  the  Uterature 
[6.  7]. 


Figure  2:  Illustration  of  the  Bowyer/Watson  algorithm  for  constructing  Delaunay  trian- 
gulations. 


In  order  to  construct  the  Delaunay  triangulation  of  a  set  of  points  using  this  algorithm,  an  initial 
triangulation  is  constructed,  either  by  forming  the  triangulation  of  the  convex  hull  of  the  points 
[2],  or  more  simply  by  creating  a  large  initial  triangulation,  using  three  or  four  auxiliary  points  (in 
two  dimensions),  which  contains  all  of  the  points  to  be  triangulated.  The  points  of  the  set  to  be 
triangulated  are  then  put  in  a  list,  and  inserted  sequentially  into  the  evolving  triangulation. 

The  Bowyer- Watson  algorithm  has  proven  to  be  very  useful  for  unstructed  mesh  generation.  It 
has  been  shunned  by  the  computational  geometry  field,  mainly  due  to  its  poor  worst-case  corirplex- 
ity  of  0(A2)  (imagine  the  case  where  each  newly  inserted  point  intersects  aU  existing  triangles). 
However,  for  unstructed  mesh  generation,  near  linear  0(N]  performance  has  been  reported  for  both 
two-  and  three-dimensional  applications  [7,  8,  9].  More  recently,  it  has  been  shown  that  the  poor 
worst-case  complexity  of  this  algorithm  represents  a  pathological  case  which  can  easily  be  avoided 
by  randomizing  the  order  in  which  the  points  are  inserted  [10]. 

2.1.2  Green-Sibson  Algorithm 

This  algorithm  is  similar  to  tho  Bowyer- Watson  algorithm  in  that  it  is  based  on  sequential  point- 
insertion  into  an  existing  triangulation,  and  also  relies  on  the  empty  circumcirde  property  [11].  To 
insert  a  new  point  into  the  triangulation,  the  triangle  which  encloses  this  new  point  is  first  located. 
The  point  is  then  inserted  into  the  triangulation  simply  by  joining  it  to  the  three  vertices  of  the 
enclosing  triangle,  as  shown  in  Figure  3.  (In  the  case  where  the  point  falls  on  a  mesh  edge,  the 
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edge  is  split  and  the  point  is  joined  to  four  mesh  vertices;  situations  in  three  dimensions  involving 
split  faces  and  edges  are  treated  analogously). 


Figure  3;  Illustration  of  the  Green-Sibson  algorithm  for  constructing  Delaunay  triangula¬ 
tions. 

The  resulting  triangulation,  although  valid,  is  not  necessarily  Delaunay,  and  the  remaining  task 
is  to  transform  it  into  the  Delaunay  triangulation  by  rearranging  the  mesh  connectivity  in  the 
vicinity  of  the  new  point.  This  is  accomplished  by  examining  the  three  (or  four)  newly  formed 
triangles.  If  their  circumcircles  are  all  empty,  then  the  triangulation  is  indeed  Delaunay,  and  no 
further  modifications  are  required.  In  the  event  one  of  these  circumcircles  contains  a  vertex,  the 
edge  of  the  triangle  which  borders  on  an  ^‘outer”  neighbor  (i.e.,  one  which  does  not  touch  the 
newly  inserted  point),  is  reconfigured  or  swapped  as  shown  in  Figure  3.  The  reason  only  these 
outer  edges  need  be  considered  is  that  the  initial  three  (or  four)  edges  which  touch  the  newly 
inserted  point  can  always  be  shown  to  be  part  of  the  final  Delaunay  triangulation  (z.e.,  consider 
how  they  may  appear  in  the  Bowyer- Watson  algorithm).  Each  time  an  edge  is  swapped,  two 
triangles  are  altered,  and  these  must  therefore  be  checked  for  the  Delaunay  criterion.  However, 
each  edge  that  is  swapped  can  be  shown  to  be  included  in  the  final  Delaunay  construction,  and 
thus  the  only  edges  which  need  be  considered  for  swapping  are  those  which  border  on  a  modified 
triangle  and  an  “outer”  untouched  triangle.  The  edge  swapping  procedure  begins  with  the  outer 
edges  of  the  newly  formed  triangles  (which  contain  the  new  point  as  a  vertex),  and  propagates 
outw^ards,  never  reexamining  the  previously  swapped  edges,  until  the  procedure  terminates  when 
no  further  edges  need  be  swapped.  There  are  certain  practical  advantages  associated  with  the 
Green-Sibson  algorithm  over  the  Bowyer- Watson  algorithm,  for  constrained  triangulations,  as  weU 
as  for  generating  tri angulations  other  than  Delaunay,  as  will  be  demonstrated. 

2.1.3  Tanemura-Merriam  Algorithm 

While  the  previous  algorithms  represent  a  top-down  approach,  a  bottom-up  approach  to  construct¬ 
ing  the  Delaunay  triangulation  of  a  given  set  of  points  is  afforded  by  the  advancing-front  Delaunay 
algorithm.  This  technique  has  apparently  been  rediscovered  several  times  in  various  fields  through¬ 
out  the  literature  [12,  13,  14,  15].  The  idea  is  to  construct  the  triangulation  one  triangle  at  a 
time,  beginning  at  the  boundaries  of  the  domain,  thus  advancing  or  sweeping  a  front  throughout 
the  domain.  The  initial  front  is  composed  of  the  set  of  edge  segments  which  define  the  convex 
hull  of  the  point-set  to  be  triangulated  (or,  for  mesh  generation  purposes,  the  set  of  edges  which 
define  aU  physical  boundaries  of  the  domain).  We  begin  by  choosing  an  edge  of  this  front.  The 
problem  consists  of  determining  the  particular  point  to  which  the  two  end  points  of  this  edge  must 
be  joined  in  order  to  construct  the  unique  Delaunay  triangle  for  this  edge,  which  will  be  present 
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ill  the  final  triaiigulation.  This  may  be  solved  in  a  simple  iterative  manner.  An  arbitrary  interior 
point  (interior  with  respect  to  the  orientation  of  the  front )  is  chosen,  and  the  triangle  formed  by  the 
two  end-points  of  the  front  edge  and  the  interior  point  is  constructed.  If  this  triangle  contains  any 
other  points  within  its  circumcircle,  it  cannot  be  a  valid  Delaunay  triangle,  and  thus  an  alternate 
point  is  chosen:  i.e.  the  point  contained  inside  the  newly  formed  circumcircle  which  is  closest  to 
its  circumcenter. 


Figure  4:  Iterative  procedure  for  determining  appropriate  point  for  constructing  the  next 
triangle  in  the  advancing-front  Delaunay  triangulation  algorithm. 


By  iterating  on  this  procedure,  as  shown  in  Figure  4,  the  appropriate  point  which  produces  a 
triangle  containing  no  points  interior  to  its  circumcircle  is  eventually  found.  This  new  triangle  is 
therefore  accepted,  and  the  front  is  advanced  by  removing  the  current  edge  from  the  front,  which 
is  now  obscured  by  the  new  triangle,  and  adding  the  new  edges  to  the  front.  (There  may  be  0,  1,  or 
2  new  edges  depending  on  whether  the  chosen  point  is  interior  to,  or  on  the  front).  The  algorithm 
terminates  once  all  edges  have  been  deleted  from  the  front,  i.e.,  when  the  entire  domain  has  been 
swept  out. 

2.1.4  Constrained  Delaunay  Triangulations 

A  triangulation  of  a  given  set  of  points  which  is  forced  to  include  as  a  subset  of  predetermined  edges, 
is  known  as  a  constrained  triangulation.  Loosely  speaking,  a  constrained  triangulation  which  is  as 
close  as  possible  to  a  Delaunay  triangulation  is  called  a  constrained  Delaunay  triangulation.  A  more 
formal  definition  of  a  constrained  Delaunay  triangulation  is  given  by  Chew  [16]:  it  is  a  triangulation 
which  contains  a  set  of  prescribed  edges,  and  such  that  the  circumcircle  of  each  triangle  contains 
no  other  vertex  of  the  mesh  which  is  visible  to  it.  A  vertex  is  visible  to  a  triangle  if  the  line  joining 
the  vertex  to  any  point  interior  to  the  triangle  does  not  intersect  one  of  the  prescribed  edges  of  the 
mesh.  An  example  of  a  constrained  Delaunay  triangulation  is  given  in  Figure  5.  The  existence  of 
constrained  Delaunay  triangulations  ensures  the  validity  of  the  Tanemura-Merriam  algorithm  for 
arbitrary  initial  (possibly  non-convex)  fronts.  It  also  guarantees  the  possibility  of  modifying  an 
existing  Delaunay  triangulation  to  include  a  set  of  edges  which  define  the  boundaries  of  the  domain 
to  be  triangulated. 
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Figure  5:  Illustration  of  constrained  Delaunay  triangulation. 

This  is  an  important  consideration  for  practical  mesh  generation  algorithms,  where  non- convex 
and  multiply  connected  domains  are  often  considered.  Unfortunately,  an  equivalent  definition  for 
constrained  Delaunay  triangulations  in  three  dimensions  is  not  available.  Thus,  the  construction 
and  modification  of  three-dimensional  Delaunay  triangulations,  which  conform  to  a  prescribed 
surface  triangulation,  has  proved  to  be  a  non-trivial  problem. 

2.2  Other  Triangulations  and  Transformation  Algorithms 

The  Delaunay  Triangulation  is  but  one  of  a  large  number  of  possible  triangulations  of  a  given  set 
of  points.  Other  possible  triangulations  include  the  minimum  total  edge-length  triangulation,  or 
minimum  weight  triangulation,  and  the  minimum  maximum-angle  triangulation  [2].  In  fact,  any 
criterion  for  describing  an  optimum  triangulation  may  be  constructed,  even  one  based  on  possible 
variable  data- values  stored  at  the  vertices  of  the  triangulation,  such  as  in  the  data-dependent  tri- 
angulations  of  [17],  and  an  algorithm  for  transforming  the  current  triangulation  to  the  so-defined 
optimum  triangulation  devised.  In  practice,  an  algorithm  capable  of  constructing  the  global  opti¬ 
mum  may  be  difficult  to  formulate,  and  local  optima,  or  improvements  from  an  initial  mesh  are 
often  produced. 

2.2.1  Two-Dimensional  Edge-Swapping 

A  method  for  transforming  a  2D  triangulation  into  another  2D  triangulation  is  given  by  the  edge- 
swapping  procedure  of  Lawson  [18].  The  algorithm  is  based  on  the  fact  that  there  are  at  most  two 
ways  of  triangulating  a  set  of  four  points,  as  shown  in  Figure  6.  For  each  pair  of  triangles  in  tne 
mesh  which  forms  a  convex  quadrilateral,  the  original  triangulation  is  compared  to  the  alternative 
triangulation  obtained  by  swapping  the  position  of  the  internal  diagonal,  as  shown  in  Figure  6.  If  the 
alternative  triangulation  is  found  to  better  optimize  the  triangulation  criterion,  then  the  diagonal  is 
swapped.  By  iterative  application  of  this  simple  diagonal  swapping  primitive  over  the  entire  mesh, 
the  triangulation  is  eventually  transformed  into  a  more  optimal  triangulation,  as  defined  by  the 
criterion.  In  [19]  it  is  proven  that  any  given  2D  triangulation  may  be  recovered  from  any  initial  2D 
triangulation  through  repeated  application  of  the  diagonal-swapping  primitive,  although  a  general 
algorithm  for  achieving  this  is  not  given.  In  general,  it  is  found  that  a  straight-forward  apphcation 
of  the  edge-swapping  procedure  gets  stuck  in  local  optima,  and  is  often  incapable  of  producing  the 
global  optimum  mesh. 


Figure  6:  Two  possible  configurations  for  the  diagonal  of  a  convex  pair  of  triangles  in  the 
edge-swapping  algorithm. 

Furthermore,  depending  on  the  order  in  which  the  edges  are  swapped,  different  local  optima  may 
be  achieved.  A  notable  exception  is  the  construction  of  a  max-min  or  Delaunay  triangulation.  It 
can  be  shown  that  repeated  apphcation  of  the  edge-swapping  technique  based  on  maximizing  the 
minimum  angle  always  converges  to  the  globally  optimum  mesh,  which  is  the  Delaunay  triangula¬ 
tion,  regardless  of  the  initial  mesh,  and  the  order  in  which  the  edges  are  swapped.  The  complexity 
of  this  procedure  is  O(NlogN),  where  N  represents  the  number  of  vertices. 

The  globally  optimum  min-max  triangulation,  on  the  other  hand,  is  not  generally  attainable 
with  this  algorithm.  A  more  sophisticated  technique,  known  as  the  edge-cutting  algorithm  [20],  is 
capable  of  transforming  any  triangulation  into  the  globally  optimum  min-max  triangulation.  The 
algorithm  may  be  interpreted  as  a  generalization  of  the  edge-swapping  procedure.  The  first  step 
consists  of  locating  the  maximum  angle  in  the  current  triangulation.  Let  prq  be  the  triangle  which 
contains  this  maximum  angle,  which  is  situated  at  point  p,  as  depicted  in  Figure  t . 


Figure  7:  Illustration  of  the  edge-cutting  algorithm  for  constructing  min-max  triangula- 
tioiis. 

In  an  attempt  to  reduce  this  angle,  a  new  point  s  is  chosen,  and  the  edge  p.s  is  drawn.  The  first 
candidate  for  the  point  s  is  the  point  of  the  neighboring  triangle  rsq.  If  the  two  angles  at  q  and  r 
are  both  smaller  than  the  original  angle  at  p,  then  the  swap  is  successful,  and  the  procedure  may 
be  continued  by  searching  for  the  next  largest  angle  in  the  mesh.  However,  if  the  angles  at  q  and  r 
are  both  larger  than  the  original  angle  at  p.  no  improvements  are  possible.  In  the  event  the  angle 
at  q  is  larger,  but  the  angle  at  r  is  smaller  than  the  original  angle  at  p,  it  can  be  shown  that  the 
optimum  new  edge  must  intersect  the  edge  rs.  Thus  we  choose  a  new  point  s’  in  the  neighboring 
triangle  of  rsq.  This  procedure  may  be  applied  recursively  to  the  neighbors  of  these  neighbors,  until 
a  suitable  point  is  found,  perhaps  several  neighbor  distances  away.  The  edge  p.s'  is  then  inserted 
into  the  triangulation,  and  all  edges  which  intersect  ps'  are  removed,  as  shown  in  Figure  i .  This 
defines  two  empty  polygons,  which  are  then  retriangulated.  This  is  achieved  by  iteratively  defining 
consecutive  vertices  of  each  polygon  boundary  which  form  a  triangle  with  all  angles  smaller  than 
the  maximum  angle  of  the  original  triangulation. 
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The  insertion  of  edges  such  as  ps',  which  may  intersect  many  existing  edges,  involves  addi¬ 
tional  non-local  information  which  enables  the  algorithm  to  avoid  getting  stuck  in  local  optima. 
However,  the  complexity  of  the  edge-cutting  algorithm  is  0{N‘^logN),  and  no  three-dimensional 
extensions  have  been  reported.  These  are  the  probable  reasons  why  this  algorithm  has  seldom  been 
implemented  for  mesh  generation  applications. 


2.2.2  Three-Dimensional  Edge-Face  S\vapping 

Three-dimensional  tetrahedralization  algorithms  based  on  simple  edge-face  swapping  primitives  are 
also  possible.  These  primitives  are  based  on  the  fact  that,  in  d  dimensions,  d  +  2  points  may  be 
triangulated  in  at  most  two  ways,  as  stated  by  Lawson  [21].  Thus,  in  three  dimensions,  a  set 
of  five  points  may  be  triangulated  in  at  most  two  manners.  Depending  on  the  configuration  of 
these  five  points,  either  a  unique  triangulation  exists,  in  which  case  the  configuration  is  called 
non-swappable,  or  two  different  triangulations  are  possible,  in  which  case  swapping  between  the  set 
of  two  constructs  is  possible.  In  general,  the  triangulation  of  five  points  in  three  dimensions  may 
result  in  2,  3,  or  4  tetrahedra,  as  shown  in  Figure  8. 


d 


d 


Figure  8:  Possible  tetrahedralizations  of  a  set  of  5  points  in  three  dimensions  (neglecting 
degenerate  cases. 


This  is  in  contrast  to  the  two-dimensional  case,  where  the  number  of  triangles  and  edges  are  identical 
for  all  possible  triangulations  of  a  given  point-set.  However,  the  triangulation  of  three-dimensional 
point-sets  no  longer  represents  a  planar  graph,  and  Euler’s  formula  no  longer  relates  the  number 
of  cells  to  the  number  of  vertices  [22].  Thus,  various  triangulations  of  the  same  point-set  in  3D 
can  be  expected  to  contain  different  numbers  of  cells,  faces,  and  edges,  and  the  face-edge  swapping 
primitives  can  be  expected  to  modify  these  numbers. 

In  the  case  where  five  vertices  are  triangulated  with  four  tetrahedra,  (i.e.,  case  (a)  of  Figure  8), 
one  of  the  five  points  is  interior  to  the  convex  hull  of  the  other  four  points,  and  no  other  triangulation 
of  this  configuration  is  possible.  This  is  a  non-swappable  configuration.  In  the  case  where  the  five 
vertices  are  triangulated  with  two  tetrahedra,  {i.e.,  case  (b)  of  Figure  8),  the  configuration  is 
swappable  provided  it  is  convex.  A  sufficient  condition  for  convexity  is  that  the  line  de  intersect 
the  face  abc  in  Figure  8  (b).  If  this  is  the  case,  then  the  two-tetrahedron  configuration  of  Figure  8 
(b)  may  be  swapped  into  the  three-tetrahedron  configuration  of  Figure  8  (c)  by  removing  face  abc 
and  inserting  the  edge  de  (and  thus  faces  ade,  bde,  cde).  If  the  initial  two-tetrahedron  configuration 
is  non-convex,  then  it  is  non-swappable.  (In  this  case,  it  corresponds  to  the  configuration  in  Figure 
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8  (c),  where  one  of  the  tetrahedra  is  missing).  Similarly,  in  the  case  where  the  five  vertices  are 
triangulated  with  three  tetrahedra,  as  in  Figure  8  (c).  the  configuration  is  swappable  into  the 
two- tetrahedron  configuration  of  Figure  8  (b),  provided  it  is  convex.  To  determine  if  the  three- 
tetrahedron  configuration  is  convex,  the  convexity  test  may  be  applied  directly  to  the  anticipated 
two-tetrahedron  swapped  configuration,  which  yields  the  same  result.  A  complete  algorithm  for 
transforming  three-dimensional  triangulations  may  be  given  as  follows; 


Step  1:  Tag  all  mesh  faces  as  candidates  for  swapping. 

Step  2:  Choose  a  pair  of  neighboring  tetrahedra  in  the  mesh  which  share  a  candidate  face.  This 
defines  five  vertices,  three  of  which  form  the  face  common  to  both  tetrahedra.  and  two  others  which 
form  the  end  points  of  the  two  tetrahedra. 

Step  3:  By  searching  through  neighboring  tetrahedra.  locate  all  other  elements  which  contain  four 
of  the  five  vertices. 

Step  4:  If  the  number  of  elements  is  four,  the  configuration  is  non-swappable.  Go  to  step  8. 

Step  5:  Check  for  convexity.  If  the  configuration  is  non-convex,  go  to  step  8. 

Step  6:  Compare  the  swapped  configuration  with  the  original  configuration.  Choose  the  configu¬ 
ration  which  “improves”  the  triangulation  based  on  the  given  criterion. 

Step  7:  Tag  the  swapped  faces  as  optimal,  and  the  neighboring  faces  as  future  candidates  for 
swapping. 

Step  8:  If  future  candidates  for  swapping  exist,  go  to  step  2,  else  end. 

As  in  the  two-dimensional  case,  these  local  transformation  techniques  most  often  terminate  in 
local  optima,  with  the  result  that  the  global  optimum  triangulation  for  a  given  criterion  is  usually 
not  achieved.  While  convergence  to  the  globally  optimum  max-min  or  Delaunay  triangulation  is 
assured  in  two  dimensions,  a  similar  result  does  not  hold  in  three  dimensions.  In  fact,  the  max- 
min  and  Delaunay  triangulation  are  not  equivalent  in  three  dimensions,  and  three-dimensional 
Delaunay  triangulations  must  be  characterized  by  the  empty-circumsphere  property.  Using  the 
empty-circumsphere  test  as  a  measure  of  optimality,  Joe  [2-3]  has  shown  that  local  transformations 
of  an  arbitrary  triangulation  are  not  guaranteed  to  converge  to  the  Delaunay  triangulations.  Only 
under  special  conditions,  such  as  the  addition  of  a  new  vertex  to  an  existing  Delaunay  triangixlation 
following  the  Green-Sibson  algorithm,  can  the  application  of  local  transformations  be  guaranteed 
to  converge  to  the  Delaunay  triangulation  [23,  24.  2.5]. 

3  Practical  Mesh  Generation  Algorithms 

The  algorithms  described  in  the  previous  section  can  be  used  to  triangulate  an  existing  set  of 
points,  or  to  modify  an  existing  triangulation.  As  such,  they  do  not  constitute  mesh  generation 
techniques,  but  rather  form  basic  building  blocks,  which  can  be  used  in  conjunction  with  other 
techniques  to  devise  mesh  generation  strategies.  In  general,  the  generation  of  unstructured  meshes 
involves  the  creation  of  both  the  mesh  points  and  their  connectivity.  This  can  be  performed  either 
sequentially  (pre-generation  of  mesh  points,  followed  by  a  triangulation  phase),  or  simultaneously 
using  an  advancing- front  technique,  or  a  Steiner  triangulation  technique.  (Steiner  triangulations 
refer  to  the  insertion  of  additional  points  into  an  existing  triangulation  in  order  to  improve  the 
quality  of  the  triangulation.)  The  problem  of  mesh  generation  includes  several  well  defined  phases 
which  can  be  summarized  as  follows: 

•  Definition  of  boundaries  of  physical  domain. 

•  Definition  of  element-size  distribution  as  a  function  of  spatial  location  (element  shape  distri¬ 
bution  is  also  required  in  the  case  of  stretched  meshes). 
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•  Generation  of  boundary  conforming  mesh  using  a  suitable  approach. 

•  Optional  mesh  post-processing  to  improve  element  quality. 

The  boundaries  of  the  physical  domain  are  usually  defined  through  some  CAD-type  data-base, 
which  most  often  involves  the  use  of  piecewise  spline  curves  in  2D,  and  assemblies  of  trimmed 
{i.e.,  possibly  intersecting)  spline  surface  patches  in  3D.  This  initial  boundary  description  must  be 
discretized  as  a  set  of  line  segments  in  2D  or  a  collection  of  planar  faces  in  3D,  for  mesh  generation 
purposes.  This  may  be  achieved  prior  to,  or  simultaneously  with,  the  construction  of  the  mesh.  The 
mesh  generation  procedure  must  be  capable  of  guaranteeing  boundary  integrity:  in  two  dimensions 
this  corresponds  to  generating  triangulations  which  contain  the  subset  of  edges  which  define  the 
domain  boundaries  (c.f.  Figure  9), 


Figure  9:  Illustration  of  non-boundary  conforming  Delaunay  triangulation  for  a  simple 
airfoil  geometry. 


while  in  three  dimensions  this  corresponds  to  generating  tetrahedralizations  which  contain  the 
subset  of  triangular  faces  which  define  the  discretized  boundary  surfaces.  The  problem  of  matching 
a  prescribed  surface  triangulation  in  the  construction  of  a  three-dimensional  volume  grid  can  be 
considerably  involved.  A  slightly  less  stringent  approach  to  three  dimensional  boundary  integrity  is 
that  of  generating  a  tetrahedralization  which  contains  a  possibly  arbitrary  boundary  triangulation 
which  nevertheless  does  not  violate  the  boundary  surface  integrity  i.e.  a  surface  triangulation  which 
does  not  intersect  or  cut  through  the  original  boundary  surface  definition.  Boundary  integrity  is 
extremely  important  in  mesh  generation,  for  if  the  boundary  surfaces  cannot  be  recovered  in  the 
mesh,  no  numerical  simulation  is  possible. 

The  definition  of  the  elemeirt-size  distribution  may  be  accomplished  either  implicitly  or  ex¬ 
plicitly.  An  implicit  definition  may  involve  inferring  a  value  for  the  element-size  in  the  interi  m 
of  the  mesh  from  the  boundary  discretization  (i.e.,  by  interpolation  from  the  closest  boundaries). 
Another  possibility  is  the  use  of  an  analytic  mapping  function  for  generating  sets  of  points  with 
desired  spacings  (i.e.,  the  use  of  structured  0-  or  C-meshes  for  generating  point  sets)  which  are 
subsequently  triangulated.  Explicit  definitions  of  the  element-size  distribution  involve  the  con¬ 
struction  of  a  function  s  =  f(x,y.z)  valid  over  the  entire  physical  domain.  This  function  may  be 
constructed  analytically,  or  using  a  small  set  of  discrete  sources  such  as: 


N 

fc=l 


1  + 


r(x)  -  Tk 


I'k 


7fc1 


(1) 


where  the  summation  is  over  all  N  sources,  and  r(x)  denotes  the  position  vector.  The  source 
parameters  Sk  and  represent  the  prescribed  element  size  at  the  source,  and  the  location  of  the 
source,  and  Tk  as  well  as  define  the  relative  region  of  influence  of  the  source. 
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Alternatively,  a  background  grid  may  be  employed,  where  at  each  vertex  a  value  of  the  ele¬ 
ment  size  is  defined,  and  values  in  between  vertices  are  obtained  by  linear  interpolation.  These 
background  meshes  may  be  unstructured  [26.  27]  (triangular,  tetrahedral),  cartesian  [28],  or  even 
quadtree-based  [29].  In  practice,  an  effective  approach  is  to  construct  an  (elliptic)  partial  differen¬ 
tial  equation  with  the  aid  of  source  terms  defined  in  a  similar  manner  to  those  of  equation  (1),  and 
to  solve  this  equation  on  the  supporting  background  grid.  The  source  terms  can  be  prescribed  in¬ 
teractively,  and  correspond  to  prescribed  mesh  spacings  at  these  physical  locations.  This  approach 
is  particularly  attractive,  since  the  elliptic  construction  ensures  a  smooth  element-size  distribution 
function,  which  is  desirable  for  the  generation  of  high  quality  meshes. 

While  mesh  quality  is  usually  defined  somewhat  heuristically,  this  usually  entails  the  notion  of 
well  shaped  elements  with  small  changes  of  element  size  between  neighboring  elements.  Regardless 
of  the  triangulation  scheme  employed  (Delaunay,  min-max,  etc..,)  it  is  now  well  known  that  accept¬ 
able  mesh  quality  cannot  be  achieved  by  optimizing  mesh  connectivity  alone,  but  is  necessarily  the 
result  of  a  close  matching  between  the  mesh-point  distribution  and  the  mesh  connectivity  employed 
in  the  construction  of  the  mesh.  (This  is  even  more  important  in  the  case  of  stretched  or  anisotropic 
meshes).  Post-processing  techniques  may  be  employed  for  improving  mesh  quality.  These  often  in¬ 
clude  minor  redistribution  of  the  mesh  points  through  Laplacian-type  smoothing  techniques,  often 
followed  by  reconfiguration  of  the  mesh  connectivity.  While  these  techniques  are  often  successful  in 
improving  grid  quality  by  removing  minor  irregularities,  they  seldom  adequately  correct  for  serious 
deficiencies  in  grid  quality.  Thus  the  ability  to  generate  a  high  quality  mesh  from  the  outset  is  of 
the  utmost  importance.  Finally,  mesh  generation  procedures  must  be  robust.  This  not  only  means 
that  they  must  use  proven  constructs  and  exact  algorithms,  but  they  must  be  designed  to  avoid 
and/or  deal  with  ambiguous  situations  which  typically  arise  in  computational  geometry  problems. 
The  use  of  exact  arithmetic  [30]  is  perhaps  the  best  approach  to  resolving  such  issues. 

3.1  Quad/Octree  Based  Mesh  Generation 

One  of  the  earliest  and  simplest  methods  for  generating  unstructured  meshes  involves  the  use  of 
quad  and  octrees  in  two  and  three  dimensions  respectively.  Considering  a  two-dimensional  example 
for  simplicity,  an  initial  quad  is  formed  which  is  large  enough  to  cover  the  entire  domain.  Assuming 
a  mesh  element-size  distribution  function  exists,  the  quadtree  is  recursively  subdivided  until  all 
leaf  quads  are  no  larger  than  the  local  value  of  the  element-size  distribution  function.  If  this 
function  is  only  defined  on  the  domain  boundaries,  the  quadtree  may  be  initially  subdivided  along 
the  boundaries.  Subdivision  in  the  domain  interior  can  then  proceed  by  ensuring  that  jumps  in 
the  sizes  of  neighboring  leaf  elements  never  exceeds  2:1.  Triangular  elements  can  be  generated 
by  forced  subdivision  of  the  quad  elements  [31],  or  by  using  the  quadtree  vertices  in  a  Delaunay 
point-insertion  algorithm  [32].  At  the  physical  boundaries,  the  quadtree  must  be  made  boundary 
conforming.  This  is  usually  accomplished  by  warping  the  mesh,  i,e,,  displacing  the  closest  quadtree 
vertices  to  coincide  with  the  boundary  curve.  Quad/octree  methods  are  relatively  simple  and 
efficient.  Their  main  deficiencies  relate  to  the  quality  of  the  mesh  near  boundaries  due  to  the 
warping  procedure. 

3.2  Advancing-Front  Methods 

Advancing- front  methods  involve  the  simultaneous  generation  of  mesh  points  and  their  connectivity. 
The  idea  is  to  build  the  mesh  element  by  element,  adding  new  elements  to  previously  generated 
elements,  thus  sweeping  out  a  front  across  the  entire  domain.  They  usually  rely  on  an  explicitly 
defined  element-size  distribution  function,  which  is  most  often  constructed  using  a  background  grid 
[26,  27,  33]. 
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Advancing- front  techniques  begin  with  a  discretization  of  the  domain  boundaries  as  a  set  of  edges 
in  two  dimensions.  This  is  accomplished  by  placing  points  along  the  boundary  curves  governed 
by  the  local  values  of  the  element-size  distribution  function,  as  shown  in  Figure  10.  These  edges 
form  the  initial  front  which  is  to  be  advanced  out  into  the  field.  A  particular  edge  of  this  front  is 
selected,  and  a  new  triangle  is  formed  with  this  edge  as  its  base,  by  joining  the  two  ends  of  the 
current  edge  either  to  a  newly  created  point,  or  to  an  existing  point  on  the  front.  The  current 
edge  is  then  removed  from  the  front,  since  it  is  now  obscured  by  the  new  triangle.  Similarly,  the 
remaining  two  edges  of  the  new  triangle  are  either  assigned  to  the  front  or  removed  from  the  front, 
depending  on  their  visibility,  as  shown  in  Figure  11. 


Figure  10:  Initial  boundary  discretization  for  the  advancing-front  method. 


The  front  thus  constitutes  a  stack  (or  priority  queue),  and  edges  are  continuously  added  to.  or 
removed  from,  the  stack.  The  process  terminates  when  the  stack  is  empty,  i.e.,  when  all  fronts 
have  merged  upon  each  other  and  the  domain  is  entirely  covered. 

The  selection  of  the  next  edge  in  the  front  may  be  based  on  various  criteria.  A  good  strategy 
is  to  always  choose  the  smallest  edge  in  the  front,  thus  ensuring  the  front  grows  from  regions  of 
small  cells  towards  regions  of  large  cells.  This  has  been  found  to  yield  smooth  high-quality  element 
distributions.  Such  a  strategy  can  easily  be  implemented  by  encoding  the  front  as  a  heap-list  of 
edges,  with  the  edge- length  as  the  heap  ordering  key  [3]. 


Figure  11:  Generation  of  new  triangle  using  new  point  (  a),  or  existing  front  point  (b),  in 
advancing-front  method. 


One  of  the  critical  features  of  such  methods  is  the  placement  of  new  points.  Upon  generating  a 
new  triangle,  a  new  point  is  first  placed  at  a  position  which  is  determined  to  result  in  an  optimal 
size  and  shape  triangle  i.e.,  along  the  median  of  the  front  edge,  at  a  distance  defined  by  the  local 
value  of  the  element-size  distribution  function,  as  shown  in  Figure  12.  The  triangle  generated  with 
this  new  point  may  result  in  a  cross-over  with  other  front  edges,  and  thus  may  be  rejected.  This  is 
determined  by  computing  possible  intersections  with  ’’nearby”  front  edges.  Alternately,  an  existing 
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point  on  the  front  may  coincidentally  be  located  very  close  to  the  new  pointy  and  thus  should  be 
employed  as  the  forming  point  for  the  new  triangle,  to  avoid  the  appearance  of  a  triangle  with  a 
very  small  edge  at  some  later  stage.  Thus,  all  front  points  within  a  certain  distance  of  the  new 
point  must  be  located.  This  involves  the  determination  of  an  appropriate  length  scale  for  defining 
the  search  region.  A  typical  example  of  a  length  scale  and  search  region  often  employed  is  shown 
in  Figure  12.  The  final  accepted  triangle  will  be  the  triangle  formed  from  the  current  front  edge 
and  the  new  point,  or  one  of  the  “nearby"  front  points,  which  does  not  intersect  any  front  edges, 
and  best  conforms  to  the  local  element-size  distribution  function. 


Figure  12:  Placement  of  new  point  Figure  13:  Possible  failure  of 

in  advancing-front  method.  advancing-front  method  for  merging 

fronts  of  dissimilar  sizes. 

The  determination  of  the  set  of  “nearby”  front  points  involves  a  proximity  search  which  is 
usually  implemented  using  a  quadtree  data-structure  [4].  Thus,  a  quadtree  based  on  the  front 
points  must  be  maintained  dynamically,  with  points  inserted  and  deleted  as  the  front  advances. 
The  front  is  thus  simultaneously  represented  as  a  heap  and  a  quadtree. 

The  space  requirements  for  such  an  algorithm  are  lower  than  may  be  expected.  Since  this 
is  essentially  a  greedy  triangulation  [2],  i,€,,  formed  elements  are  never  subsequently  modified, 
all  points,  edges  and  triangles  which  lie  behind  the  front  need  no  longer  be  considered  in  the 
generation  process.  Thus  the  only  active  portion  of  the  data  is  the  front.  Since  a  front  has  one 
lower  dimension  than  the  domain  to  be  discretized,  the  required  space  for  such  an  algorithm  in 
two  dimensions  is  0(\/iV),  where  N  is  the  final  number  of  grid  points  generated.  Since  N  points 
are  added  sequentially,  the  complexity  is  at  most  0{N\/N).  However,  by  employing  sophisticated 
searching  techniques  such  as  spatial  quad-trees,  this  complexity  is  easily  lowered  to  0{Nlog\/N) 
v/hicli  is  asymptotically  equivalent  to  0{NlogN).  Optimal  space  usage  has  not  in  general  been 
achieved,  due  to  the  difficulty  in  continuously  dumping  out  generated  elements.  However,  restart 
capabihties  are  easily  implemented  [34,  29].  which  can  greatly  reduce  the  required  working  size  for 
a  large  mesh  generation  job. 

One  of  the  advantages  of  such  an  approach  is  the  automatic  point  placement  strategy,  which 
generally  results  in  high-quality  elements  throughout  most  of  the  flow-field,  due  to  the  optimum 
positioning  of  these  new  points.  Additionally,  all  real  operations  performed  (such  as  intersection 
checking)  are  of  a  local  nature;  f.e.,  intersection  checks  are  performed  with  neighboring  edges  of 
similar  length,  thus  reducing  the  chances  for  round-off  error  induced  failure.  Finally,  boundary 
integrity  is  guaranteed,  since  the  boundary  discretization  constitutes  the  initial  condition. 

The  disadvantages  of  advancing-front  techniques  mainly  relate  to  their  efficiency.  The  inter¬ 
section  checking  phase  is  a  rather  brute-force  technique  for  ensuring  the  acceptability  of  a  new 


triangle,  which  is  relatively  expensive.  Additionally,  for  each  generated  triangle,  the  quad-tree 
data-structure  must  be  traversed  from  top  to  bottom  (0{logN)  steps)  in  order  to  locate  ’’nearby” 
points  and  edges.  Another  contributing  factor  is  the  fact  that  advancing-front  techniques  construct 
the  mesh  one  triangle  at  a  time.  Since  in  two  dimensions  there  are  asymptotically  twice  as  many 
triangles  as  points,  a  more  efficient  strategy  would  be  to  construct  the  mesh  one  point  at  a  time. 
Thus,  each  time  a  new  point  is  generated,  efficiency  could  be  gained  by  determining  all  the  potential 
triangles  which  join  this  new  point  to  the  existing  front  with  a  single  traversal  of  the  quad-tree 
data-structure.  In  three  dimensions,  the  savings  are  even  greater,  since  there  exists  on  average  .5 
to  6  times  more  tetrahedra  than  vertices. 

Finally,  even  though  advancing-front  techniques  rely  only  on  local  operations,  they  may  still 
suffer  from  robustness  problems.  Central  to  the  issue  of  determining  acceptable  triangles  and  ’’best” 
points,  is  the  determination  of  a  local  length  scale  which  defines  the  region  of  ’’nearby”  points 
and  edges.  This  length  scale  is  generally  obtained  from  the  field  function  (which  may  employ  a 
background  grid).  In  the  case  of  two  merging  fronts,  if  the  field  function  varies  rapidly  over  the 
region  between  the  merging  fronts,  the  relative  sizes  of  the  edges  on  one  front  may  be  much  larger 
than  those  on  the  other  front.  If  a  search  is  initiated  from  the  front  with  the  smaller  length  scale, 
the  region  of ’’nearby”  edges  may  not  contain  the  appropriate  edges  and  points  of  the  other  front, 
and  failure  will  occur,  as  shown  in  Figure  13.  Thus,  the  success  of  the  advancing-front  technique 
relies  on  the  existence  of  a  smoothly  varying  field  function. 

3.3  Delaunay  Point-Insertion  Methods 

Delaunay-based  methods  offer  the  possibility  of  constructing  mesh  generation  strategies  based 
on  proven  computational-geometry  algorithms  (although  this  does  not  in  itself  guarantee  a  more 
efficient/robust  overall  approach). 

Some  of  the  earhest  Delaunay-based  mesh  generation  strategies  relied  on  predetermined  mesh- 
point  sets  [8,  7,  3.5,  36].  For  multi-component  geometries,  each  component  represents  a  simple 
configuration  which  can  be  fitted  locally  with  a  simple  structured  mesh  {i.e.,  for  example  using 
a  structured  0-mesh  about  each  element  of  a  multi-element  airfoil.  By  constructing  a  set  of 
overlapping  structured  meshes  in  this  manner,  and  discarding  the  connectivity  of  these  meshes,  as 
well  as  the  points  which  fall  outside  of  the  physical  domain,  a  set  of  points  is  obtained  which  may  be 
used  as  the  basis  for  a  Delaunay  triangulation  mesh  generation  strategy,  using  the  Bowyer/Watson 
or  Green-Sibson  algorithms. 

In  one  particular  approach  [8,  3-5],  an  initial  Delaunay  mesh  is  constructed  by  joining  one  of  the 
inner  boundary  points  to  all  of  the  outer  boundary  points  in  a  hub  and  spoke  type  arrangement. 
The  inner  boundary  points  are  then  inserted  into  the  triangulation  using  the  Bowyer/Watson  algo¬ 
rithm.  After  all  inner  boundary  points  are  inserted,  a  check  for  boundary  integrity  is  performed.  If 
boundary  violations  are  encountered,  additional  boundary  mesh-points  can  be  inserted  to  remedy 
the  situation.  Once  boundary  integrity  has  been  recovered,  triangles  exterior  to  the  domain  are 
removed,  thus  effectively  freezing  all  boundary  defining  edges,  and  all  interior  points  are  inserted 
and  triangulated  using  the  Bowyer/Watson  algorithm.  Due  to  the  overlapping  mesh  construction 
of  the  point-set,  coincidentally  close  mesh  points  may  be  produced,  and  mesh  smoothing  is  em¬ 
ployed  as  a  post-processing  operation  in  order  to  relieve  effects  caused  by  local  irregularities  in  the 
final  mesh  spacing,  as  seen  in  Figure  14.  Mesh-point  filtering  techniques  may  also  been  employed 
to  restore  smooth  mesh-point  distributions  [37].  These  techniques  have  also  been  employed  for 
three-dimensional  mesh  generation  about  aircraft-type  configurations  by  Baker  [7,  38].  Although 
predetermined  mesh-point  set  techniques  have  been  very  successful  for  certain  classes  of  geome¬ 
tries,  the  simultaneous  mesh-point  generation  and  triangulation  techniques  embodied  in  Steiner 
triangulations  offer  improved  flexibility  and  automation  for  arbitrary  geometries.  A  Steiner  tri¬ 
angulation  is  a  triangulation  whose  additional  points  are  inserted  in  order  to  improve  the  quality 
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of  the  triangulation.  Assuming  an  initial  triangulation  and  an  element-size  distribution  function 
exist,  the  triangulation  can  be  incrementally  modified  by  continually  inserting  new  points  until  the 
final  mesh  closely  matches  the  prescribed  element-size  distribution  function  in  all  regions  of  the 
domain.  A  particularly  effective  strategy  is  to  insert  the  new  mesh  points  at  the  circumcenters  of 
the  triangles  which  are  flagged  for  improvement  [39,  40,  41,  9,  25]. 


Figure  14:  Delaunay  triangulation  of  mesh-point  distribution  generated  from  overlapping 
structured  meshes,  before  and  after  application  of  smoothing. 


This  strategy  can  be  proved  to  result  in  triangulations  where  the  angles  are  aU  bounded  between 
30^  and  120'^,  neglecting  boundary  effects  [16,  42].  This  result  is  a  consequence  of  the  fact  that  the 
circumcenter,  where  the  new  point  is  placed,  is  equidistant  from  the  three  vertices  of  its  forming 
triangle,  which  represent  the  closest  points  to  the  new  point,  since  the  circle  is  necessarily  empty  by 
the  Delaunay  criterion.  An  effective  mesh  generation  strategy  using  this  approach  can  be  formulated 
as  follows: 


Step  1:  C  ■onstruct  an  element-size  distribution  function  .5  =  f{x,y^z). 

Step  2:  Discretize  all  boundary  curves  based  on  the  above  function. 

Step  3:  Construct  an  initial  triangulation  which  covers  the  entire  domain  (usually  by  forced 
triangulation  of  a  large  quadrilateral  in  2D  or  hexahedron  in  3D). 

Step  4:  Insert  all  the  boundary  points  into  the  triangulation  using  the  Bowyer/Watson  or  Green- 
Sibson  algorithm. 

Step  5:  Construct  a  heap-list  of  all  triangles  in  the  mesh  which  are  larger  than  the  local  value 
specified  by  the  element-size  distribution  function  (usually  all  existing  triangles  at  this  initial  stage); 
the  heap-list  is  ordered  by  some  measure  of  triangle  size  (f.e.,  circumradius)  or  quality  (maximum 
angle). 

Step  6:  PuU  the  first  triangle  off  the  top  of  the  heap- list,  insert  a  new  point  at  its  circumcenter  and 
retriangulate  using  the  Bowyer/Watson  or  Green-Sibson  algorithm.  For  each  of  the  newly  formed 
triangles,  which  are  larger  than  the  local  value  specified  by  the  element-size  distribution  function, 
insert  them  into  the  heap-list. 

Step  7:  If  the  heap-list  is  empty,  stop,  else  go  to  step  6. 

Step  8:  Recover  the  boundary  integrity. 

The  heap-list  empties  out  when  the  triangulation  converges  to  the  element-size  distribution  pre¬ 
scribed  by  the  function  ^  =  f{x,y,z).  This  method  results  in  a  very  efficient  mesh  generation 
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technique.  Since  the  identity  of  the  originating  triangle  for  the  circumcenter  point-insertion  opera¬ 
tion  is  initially  known,  the  search  for  the  containing  triangle  in  the  Delaunay  algorithm  is  initiated 
with  this  triangle.  In  most  cases,  this  corresponds  to  the  containing  triangle,  and  in  all  other  cases, 
the  containing  triangle  is  only  several  neighbors  away.  Thus,  the  union  of  intersected  triangles 
in  Watson’s  algorithm  can  usually  be  found  very  rapidly,  in  constant  time.  Since  the  heap-list 
operations  (insert/delete)  can  all  be  performed  in  OilogN)  time,  the  space  and  time  requirements 
of  this  algorithm  are  0{NlogN).  In  practice,  this  method  has  produced  some  of  the  most  efficient 
mesh  generation  codes  available  today,  capable  of  generating  1  million  three-dimensional  elements 
per  hour  on  present-day  workstations  [43]. 

In  the  above  description  of  the  algorithm,  the  boundary  recovery  procedure  is  performed  as  a 
final  step,  after  all  mesh  points  have  been  inserted.  In  practice,  this  procedure  may  be  performed 
either  in  the  final  stages  of  the  algorithm,  or  earlier  on,  just  after  the  boundary  points  have  been 
inserted,  i.e.,  between  steps  4  and  5.  Recovering  the  boundary  integrity  in  the  final  stages  of  the 
algorithm  enables  the  boundaries  to  be  neglected  and  the  domain  to  be  treated  as  convex,  during 
the  point-insertion  phase.  This  permits  a  straight-forward  application  of  the  Delaunay  insertion 
algorithm  at  this  stage.  Furthermore,  mesh  quality  is  much  higher  at  the  end  of  the  point-insertion 
procedure,  which  may  result  in  a  more  robust  boundary  recovery  phase.  On  the  other  hand,  new 
points  may  be  introduced  arbitrarily  close  to  the  eventual  boundary  surfaces,  during  the  point- 
insertion  process,  which  may  result  in  a  poor  quality  mesh  after  the  boundary  recovery  operation. 

If  the  boundary  recovery  is  performed  just  after  the  insertion  of  the  boundary  points,  interior 
mesh  points  near  the  boundary  may  be  positioned  more  optimally  (see  for  example  section  3.5). 
However,  this  requires  boundary  integrity  be  maintained  subsequently  throughout  the  remainder 
of  the  grid  generation  process,  and  implies  the  application  of  Delaunay  point-insertion  algorithms 
in  non-con  vex  domains  for  step  6.  In  two  dimensions,  the  existence  of  a  constrained  Delaunay 
triangulation  guarantees  that  the  point-insertion  algorithms  can  be  extended  by  terminating  the 
search  for  intersected  triangles  in  the  given  direction,  when  a  boundary  is  encountered.  In  three 
dimensions,  although  similar  techniques  have  been  employed,  there  is  no  guarantee  that  a  valid  tri¬ 
angulation  will  result,  since  no  three-dimensional  definition  of  a  constrained  Delaunay  triangulation 
is  available. 

One  way  of  recovering  the  boundary  integrity,  either  at  the  end  of  the  mesh  generation  process, 
or  after  the  insertion  and  triangulation  of  the  prescribed  boundary  points,  is  to  insert  a  sufficient 
number  of  additional  boundary  points  until  the  Delaunay  triangulation  conforms  to  the  boundary. 
A  more  elegant  approach  consists  of  modifying  the  existing  Delaunay  triangulation  using  the  edge 
swapping  (in  2D)  and  face-edge  swapping  (in  3D)  primitives  discussed  in  section  2.2  to  produce 
a  locally  modified  triangulation  which  conforms  to  the  boundary.  This  approach  also  offers  the 
possibility  of  matching  exactly  a  predetermined  boundary  discretization.  In  two  dimensions,  each 
edge  of  the  prescribed  boundary  discretization  is  searched  for,  in  the  Delaunay  mesh.  If  one  of 
these  edges  is  not  found,  but  the  two  endpoints  of  the  edge  are  located  in  neighboring  triangles, 
then  a  simple  edge  swap,  as  depicted  in  Figure  6,  is  sufficient  to  recover  the  missing  edge.  In  the 
event  these  endpoints  are  located  in  non-neighboring  triangles,  then  aU  triangles  which  the  missing 
edge  intersects  are  first  located.  These  triangles  are  aU  removed,  and  the  missing  edge  is  inserted, 
as  shown  in  Figure  15.  This  creates  two  polygons,  on  either  side  of  the  new  edge,  which  are  then 
triangulated,  using,  for  example,  the  Tanemura-Merriam  algorithm. 

The  situation  in  three  dimensions  is  somewhat  more  complex.  The  procedure  is  usually  divided 
into  two  phases  for  a  prescribed  boundary  surface  triangulation:  firstly  the  recovery  of  the  boundary 
edges,  and  secondly  the  recovery  of  the  boundary  faces.  The  face-edge  swapping  primitives  of  section 
2.2  are  utibzed  to  transform  the  Delaunay  triangulation  into  one  which  contains  the  boundary 
edges  and  faces.  However,  contrary  to  the  two-dimensional  case,  the  insertion  of  additional  points 
is  often  required  to  enable  the  completion  of  the  boundary  recovery  procedure  in  three  dimensions. 
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VVeatherill  [40]  lias  categorized  all  possible  boundary  edge  and  face  intersection  types,  and  lists 
rules  for  inserting  new  points  to  enable  recovery  of  these  boundary  elements.  Generally  speaking, 
he  advocates  the  creation  of  new  points  on  the  surface,  at  the  point  of  intersection  with  the  existing 
grid  element.  Once  all  boundary  edges  and  faces  have  been  recovered,  the  new'  surface  points  may 
be  removed  in  an  attempt  to  recover  the  original  surface  discretization.  However,  in  some  cases, 
the  removal  of  these  point  may  result  in  cavities  which  cannot  be  re-tetrahedralized. 


Figure  15:  Illustration  of  the  insertion  of  a  missing  boundary  edge  into  an  existing  De¬ 
launay  triangulation,  and  the  resulting  empty  polygons  which  must  be  retriangulated. 

In  [9],  a  similar  method  wdiich  introduces  new  points  off  the  surface  is  developed,  thus  permitting 
the  recovery  of  the  original  surface  discretization. 

The  mesh-point  distribution  produced  by  the  circumcenter  point-placement  strategies,  while 
certainly  acceptable,  are  somewhat  irregular  and  lack  the  high  degree  of  smoothness  produced 
by  the  advancing-front  method.  This  is  perhaps  due  to  the  fact  that  Steiner  triangulations  are 
top-down  approaches,  wdiich  seek  to  improve  existing  triangulations  through  refinement.  A  more 
serious  difficulty  with  Delaunay  point-insertion  methods  relates  to  robustness  problems  due  to 
round-off  error.  Particularly  in  the  initial  phases,  when  the  boundary  points  are  inserted,  highly 
distorted  triangular /tetrahedral  elements  are  formed.  The  finite  precision  real  arithmetic  utilized 
to  compute  all  quantities  [  i.e.,  cell  volume,  circumcenter,  etc.)  may  lead  to  failure  of  the  algorithm 
simply  due  to  insufficient  accuracy  (even  in  64  bit  arithmetic).  One  remedy  is  to  alter  the  order 
in  which  the  points  are  inserted,  or  to  employ  auxiliary  (possibly  temporary)  points  to  avoid  the 
creation  of  highly  skewed  elements.  How'ever,  this  may  also  have  the  effect  of  modifying  the  resulting 
triangulation.  Ultimately,  the  best  technique  for  ensuring  robustness  is  to  resort  to  infinite  precision 
arithmetic  using  integer  programming  techniques  [30],  perhaps  in  combination  with  some  of  the 
above  strategies. 

3.4  Advancing-Front  Delaunay  Triangulation 

In  an  effort  to  alleviate  some  of  the  drawbacks  of  point-insertion  Delaunay  triangulation  methods,  an 
advancing-front  Delaunay  triangulation  algorithm  was  recently  developed  by  the  author  [29].  This 
method  enjoys  the  smooth  point  distribution  and  guaranteed  boundary  integrity  of  the  traditional 
advancing-front  methods,  while  using  the  well  founded  principles  of  Delaunay  triangulation  to 
replace  the  heuristics  present  in  the  reconnection  phase  of  traditional  advancing-front  methods. 
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Robustness  without  the  use  of  high  accuracy  arithmetic  is  assured  by  only  constructing  triangular 
elements  which  are  acceptable  in  size  and  shape,  as  determined  by  the  background  spacing  function, 
i.e.,  the  highly  skewed  triangles  which  appear  temporarily  in  the  point-insertion  methods  are  never 
formed. 

A  background  spacing  function  s  =  f{x,y)  is  used  to  determine  the  maximum  permissible 
circumradius  p  of  a  triangle  as  a  function  of  spatial  location.  Each  time  a  new  point  is  added  ahead 
of  the  front,  it  is  desired  to  construct  all  Delaunay  triangles  which  contain  this  new  point,  but  which 
do  not  violate  the  local  circumradius  bound.  Triangles  which  violate  the  local  circumradius  bound 
should  not  be  constructed,  even  temporarily,  for  this  may  require  non-local  operations  and  the 
possibility  of  round-off  induced  error.  One  method  of  constructing  these  triangles  is  simply  to  join 
the  new  point  to  every  possible  pair  of  points  in  the  grid  and  preserve  each  potential  triangle  which 
does  not  violate  the  Delaunay  criterion  and  the  circumradius  bound.  A  more  efficient  technique 
is  to  determine  a  subset  of  the  grid  points  which  is  sufficient  for  locating  aU  acceptable  triangles. 
Such  a  subset  can  be  formed  by  considering  all  front  points  which  are  less  than  2p  away  from  the 
new  point.  The  advancing-front  Delaunay  triangulation  algorithm  for  fixed  point  sets  (Tanemura- 
Merriam)  described  in  section  2.1.3  [14,  15]  can  then  be  employed  to  construct  the  Delaunay 
triangles  between  these  front  points  and  the  new  point.  As  in  the  traditional  advancing-front 
algorithm,  it  is  possible  that  a  new  triangle  may  be  formed  from  the  current  front  edge  with  an 
existing  front  point.  Such  cases  are  automatically  detected  by  the  Tanemura-Merriam  algorithm, 
and  the  new  point  is  rejected.  Finally,  when  a  new  point  is  introduced,  it  is  possible  that  it  intersects 
circumcircles  of  existing  triangles.  If  this  is  the  case,  these  triangles  must  be  deleted  prior  to  the 
formation  of  new  triangles,  since  they  are  in  violation  of  the  Delaunay  criterion.  The  method  is 
set  up  similarly  to  the  traditional  advancing-front  procedure.  A  background  spacing  function  is 
defined,  and  the  boundaries  are  discretized  and  constitute  the  initial  front.  The  basic  algorithm 
can  be  summarized  as  follows: 


Step  1:  Construct  the  original  front  as  a  set  of  boundary  edges. 

Step  2:  Choose  a  particular  edge  from  the  front,  according  to  some  criteria  such  as  minimum  edge 
length. 

Step  3:  Locate  all  front  points  which  are  less  than  2p  away  from  either  end  point  of  this  edge, 
where  p  is  the  local  circumradius  bound  determined  from  the  field  function. 

Step  4:  Use  the  Tanemura-Merriam  algorithm  to  determine  the  Delaunay  triangle  formed  between 
this  edge  and  the  set  of  candidate  points,  if  such  a  triangle  exists. 

Step  5:  If  this  triangle  exists  and  is  acceptable  (circumradius  smaller  than  p),  form  a  new  triangle, 
update  the  front,  and  proceed  to  step  12.  Otherwise  create  a  new  point  at  the  appropriate  location. 
Step  6:  Determine  aU  front  triangles  whose  circumcircles  are  intersected  by  the  new  point. 

Step  7:  Using  a  neighbor  search  initiated  at  the  intersected  front  triangles,  locate  all  interior 
triangles  whose  circumcircles  are  intersected  by  the  new  point. 

Step  8:  Jlemove  all  such  triangles  and  update  the  front.  Any  new  front  points  which  result  from 
this  operation  are  added  to  the  list  of  “close”  points. 

Step  9:  If  the  circumradius  of  any  of  the  new  intersected  triangles  is  larger  than  the  previously 
determined  maximum  permissible  value  p,  replace  the  old  value  by  this  new  maximum,  and  locate 
any  additional  front  points  which  are  less  than  2p  away  from  the  new  point. 

Step  10:  Form  aU  possible  Delaunay  triangles  which  contain  the  new  point  and  two  other  points 
in  the  list  of  “close”  points,  and  which  do  not  violate  the  local  circumradius  bound.  Such  triangles 
are  found  using  the  Tanemura-Merriam  algorithm. 

Step  11:  Add  these  triangles  to  the  mesh  and  update  the  front. 

Step  12:  If  the  front  queue  is  empty,  stop,  otherwise  go  to  step  2. 


19 


As  in  previous  advanciiig-front  algorithms,  the  front  is  stored  as  a  heap-list  in  order  to  facilitate 
the  choice  in  step  2.  A  front-based  quadtree  is  also  used  to  perform  the  search  for  ‘"close”  points  in 
step  3.  A  second  quadtree  is  also  employed  in  the  search  for  intersected  front  triangle  circumcircles. 
(The  construction  of  quadtrees  for  geometric  entities  other  than  points,  such  as  circles,  is  described 
in  [4]). 

In  most  cases  when  a  new  point  is  inserted  ahead  of  the  front,  no  intersected  circumcircles  are 
detected,  and  few  if  any  other  front  points  are  within  the  critical  2p  distance  of  the  new  point,  thus 
the  construction  of  new  triangles  is  quite  simple.  When  fronts  merge  upon  each  other,  the  situation 
becomes  more  complex.  Triangle  circumcircles  from  the  opposing  front  are  usually  intersected.  This 
alerts  the  front  being  advanced  to  the  proximity  of  the  neighboring  front  and  provides  an  extra 
length  scale  (the  circumradius  of  the  intersected  triangle)  which  is  utilized  effectively  in  mending 
the  two  fronts  as  shown  in  Figure  16. 


Figure  16:  Illustration  of  the  use  of  the  front  triangle  circumcircle  to  detect  the  proximity 
of  a  front,  (compare  with  the  situation  depicted  in  Figure  13). 


The  Delaunay  construction  naturally  results  in  additional  information  which  is  essential  for  merging 
two  fronts  of  dissimilar  length  scales.  This  is  precisely  the  information  which  is  missing  in  the 
traditional  advancing-front  algorithms,  which  are  prone  to  failure  when  fronts  of  dissimilar  scales 
are  encountered.  Although  a  smoothly  varying  background  function  can  reduce  the  likelihood  of 
having  to  merge  dissimilar  fronts,  this  cannot  be  guaranteed.  The  extra  information  inherent  in 
the  Delaunay-based  algorithm  permits  the  merging  of  arbitrarily  dissimilar  fronts,  thus  increasing 
robustness. 

The  space  requirements  and  computational  efficiency  of  the  present  algorithm  lie  in  between 
those  of  traditional  advancing-front  algorithms  and  the  Bowyer/Watson  algorithm  for  Delaunay  tri¬ 
angulation.  As  opposed  to  the  advancing-front  algorithms,  the  present  approach  does  not  represent 
a  true  greedy  algorithm  [2],  he.,  triangles  behind  the  front  may  be  subsequently  modified.  However, 
the  only  such  triangles  which  may  be  modified  are  those  whose  circumcircle  extends  ahead  of  the 
front  into  the  ungridded  region  into  which  new  points  are  placed.  Assuming  a  relativ^ely  smooth 
distribution  of  elements  behind  the  front,  the  number  of  such  non-frozen  elements  is  a  constant 
times  the  size  of  the  front.  Thus,  we  can  expect  a  space  requirement  of  0{  \/N),  although  the  worst 
case  estimate  is  more  likely  0(N).  On  the  other  hand,  it  is  a  simple  matter  to  create  a  restart 
facility  which  dumps  out  the  generated  portion  of  the  grid  after  a  prescribed  number  of  elements 
have  been  produced,  and  reinitializes  the  generation  process  using  the  front  of  the  previous  mesh 
as  the  initial  condition.  If  no  old  elements  behind  the  front  are  considered  in  the  restart  process, 
the  resulting  mesh  may  contain  regions  of  locally  non-Delaunay  triangles  along  the  fronts  present 
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at  each  restart  phase.  If  a  true  Delaunay  triangulation  is  required,  these  regions  may  be  converted 
using  the  edge-swapping  algorithm  in  a  postprocessing  phase. 

The  current  algorithm  exhibits  a  worst  case  complexity  of  0{N‘^),  just  as  the  Bowyer/Watson 
algorithm  for  Delaunay  triangulation.  This  occurs  when  the  circumcircles  of  all  existing  triangles 
are  intersected  by  each  new  point,  or  when  all  front  points  must  be  included  in  the  list  of ’’nearby” 
points  which  are  candidates  for  forming  a  new  element.  However,  for  the  smooth  element  and  point 
distributions  which  are  sought  in  the  context  of  mesh  generation,  the  number  of  points  within  the 
characteristic  distance  of  a  newly  inserted  point  and  the  number  of  intersected  triangles  should 
approach  a  constant.  When  the  logN  term  from  the  quad-tree  data-structures  employed  for  the 
search  routines  on  the  front  is  included,  a  complexity  of  O(NlogN)  can  be  expected.  This  is 
the  same  complexity  exhibited  by  other  advancing-front  algorithms  under  the  same  assumptions. 
However,  the  present  algorithm  can  be  expected  to  run  significantly  faster  than  other  advancing- 
front  algorithms,  since  the  mesh  is  generated  one  point  at  a  time,  rather  than  one  triangle  at  a 
time.  In  two  dimensions,  the  differences  may  be  small,  especially  since  two  length  scales  and  thus 
two  searches  on  the  front  are  required  for  robustness  (an  additional  one  for  the  intersected  front 
triangle  circumcircles).  However,  in  three  dimensions  where  there  are  on  the  average  5  to  6  times 
more  tetrahedra  than  vertices,  the  0{logN)  cost  of  traversing  the  octree  data-structures  may  be 
amortized  over  aU  elements  generated  about  each  newly  inserted  mesh  point. 


Figure  17:  Illustration  of  advancing-front  Delaunay  triangulation  procedure,  (i)  Quadtree 
of  initial  boundary  discretization,  (ii)  Grid  at  intermediate  stage  of  generation,  (iii)  Final 
generated  mesh  prior  to  the  application  of  smoothing,  (iii)  Final  smoothed  mesh. 


On  the  other  hand,  the  present  algorithm  will  probably  not  achieve  the  efficiency  exhibited  by 
Delaunay-triangulation  point-insertion  methods,  due  to  the  need  to  traverse  the  quad-tree  data- 
structures,  which  are  not  present  in  these  other  methods,  and  the  need  to  consider  a  sufficient  but 
not  necessary  list  of  candidate  points  for  triangulation  at  each  point-insertion  process.  This  cost,  as 
well  as  the  increased  coding  complexity,  is  viewed  as  the  price  required  for  additional  robustness.  An 
example  of  the  mesh  generation  process  using  the  advancing-front  Delaunay  triangulation  algorithm 
is  depicted  in  Figure  17.  The  initial  front  quadtree  is  shown  in  the  first  part  of  the  figure.  Two  copies 
of  this  quadtree  are  kept.  The  first  is  static  and  forms  the  supporting  structure  for  the  background 
spacing  function.  In  order  to  evaluate  this  function  at  a  given  spatial  location  .s  f{x,y),  the 
quadtree  is  descended  to  the  leaf  quad  which  contains  the  particular  (x,  y)  location,  and  the  function 
value  is  obtained  by  bilinear  interpolation  of  the  values  at  the  four  quad  vertices.  The  second  copy 
of  this  quadtree  is  dynamic,  and  is  advanced  along  with  the  front  throughout  the  mesh  generation 
process.  This  quadtree  is  employed  for  the  spatial  searches  required  by  the  algorithm.  The  second 
part  of  the  figure  illustrates  an  intermediate  stage  of  the  mesh  generation  process,  while  the  final 
mesh  is  depicted  in  the  third  part  of  the  figure,  and  the  smoothed  mesh  in  the  fourth  part  of 
the  figure.  It  is  noteworthy  that  the  final  mesh  generated  in  this  fashion  exhibits  a  very  regular 
distribution  of  elements,  even  before  the  application  of  post-processing  techniques.  (Compare  this 
with  the  first  part  of  Figure  14,  for  example). 

3.5  Advancing-Front  Point-Insertion  Methods 

Another  option  for  combining  the  advantages  of  Delaunay  triangulation  with  those  of  advancing- 
front  techniques  is  to  employ  the  point-insertion  algorithms  (e.g.  Bowyer/Watson  or  Green-Sibson) 
in  conjunction  with  an  advancing-front  type  point-placement  strategy.  The  advantages  of  such 
methods  are  the  smooth  point  distributions  associated  with  advancing-front  methods,  coupled 
with  the  efficiency  and  ease  of  implementation  of  the  Delaunay  point-insertion  methods.  The 
disadvantages  are  the  same  as  those  described  for  the  circumcircle  point-insertion  methods  of 
section  3.3,  i.e.,  the  application  of  boundary  recovery  techniques,  and  possible  robustness  problems 
associated  with  the  use  of  finite-precision  arithmetic,  due  to  the  presence  of  highly  skewed  elements 
at  intermediate  stages  of  the  mesh  generation  process.  However,  these  methods  are  considerably 
simpler  to  implement  than  the  algorithm  of  the  previous  section.  They  extend  readily  to  three 
dimensions,  and  can  also  be  used  to  generate  triangulations  other  than  the  Delaunay  tri angulation, 
as  will  be  shown.  Such  methods  have  been  described  by  Rebay  [44],  and  Mueller  et  al.  [45],  and 
later  taken  up  by  Marcum  et  al.  [43,  46].  Assuming  an  element-size  distribution  function  has  been 
defined,  and  is  used  to  control  the  maximum  permissible  circumradius  of  a  triangle  as  a  function  of 
spatial  location,  Rebay  [44]  begins  by  classifying  all  triangles  of  the  initial  tri  angulation  as  either 
accepted,  active,  or  waiting.  An  accepted  triangle  is  one  that  satisfies  the  circumradius-bound 
defined  by  the  field  function,  while  a  waiting  triangle  is  one  that  does  not.  Active  triangles  are 
waiting  triangles  which  are  neighbors  of  accepted  triangles.  Active  triangles  are  the  only  triangles 
considered  for  refinement.  Thus,  the  set  of  edges  which  delimits  the  border  between  accepted  and 
waiting  triangles  constitutes  the  front,  and  active  triangles  are  simply  triangles  which  have  an 
edge  on  the  front.  Mueller  et  al,  [45]  and  Marcum  [46]  have  devised  point  placement  strategies 
which  mimic  hyperbolic  structured  mesh  generation  techniques.  At  each  stage,  the  entire  front  is 
considered,  and  a  new  set  of  points  is  created  ahead  of  the  front,  similarly  to  a  new  layer  of  points 
in  a  structured  hyperbolic  mesh  scheme.  These  points  are  then  filtered,  in  order  to  remove  points 
which  may  be  too  close  together,  such  as  in  concave  regions  of  the  front,  and  the  filtered  point-set  is 
inserted  into  the  mesh  sequentially,  using  the  Delaunay  point-insertion  methods.  A  more  traditional 
technique  for  advancing  the  front  by  using  a  heap-list,  which  always  chooses  the  smallest  front-edge, 
has  also  been  found  to  work  well,  A  particular  implementation  of  these  algorithms,  as  performed 
by  the  author  is  given  below: 
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step  1:  Define  afield-spacing  function  and  discretize  the  boundaries  accordingly. 

Step  2:  Create  an  initial  triangulation  which  covers  the  entire  domain  (using  auxiliary  points) 
Step  3:  Insert  all  boundary  points  into  the  triangulation  using  a  Delaunay  point-insertion  algo¬ 
rithm  (Bowyer/Watson  or  Green- Sibson). 

Step  4:  Recover  the  boundary  integrity  using  edge-swapping  techniques  (or  edge-face  swapping  in 
3D) 

Step  5:  Flag  all  triangles  outside  of  the  physical  domain  as  accepted,  and  all  triangles  inside  the 
domain  which  violate  the  local  circumradius-bound  as  waiting  (usually  all  initial  triangles  will  be 
in  violation). 

Step  6:  The  set  of  edges  which  delimits  accepted  and  waiting  triangles  forms  the  current  front 
(which  coincides  with  the  boundaries  at  this  stage).  Construct  a  heap-list  of  these  edges  ordered 
by  smaRest  edge-length. 

Step  7:  Pick  the  smallest  front  edge  from  the  top  of  the  heap-list,  and  determine  the  position  of 
the  new  point  ahead  of  this  edge,  based  on  the  local  value  of  the  field  function. 

Step  8:  Insert  this  point  into  the  triangulation  using  the  Delaunay  point-insertion  algorithm. 
Step  9:  Reclassify  newly  formed  triangles  as  accepted  and  waiting,  based  on  the  local  circumradius 
bound,  and  update  the  front  accordingly. 

Step  10:  If  the  front  is  not  empty,  go  to  step  7,  else  stop. 


Figure  18:  Illustration  of  point-placement  strategy  for  advancing-front  Delaunay  point- 
insertion  method. 

The  point -placement  strategy  is  illustrated  in  Figure  18.  The  new  point  is  placed  along  the  median 
of  the  current  front-edge,  at  a  distance  which  results  in  a  triangle  of  circumradius  p,  as  prescribed 
by  the  field  function.  The  location  of  the  new  point  along  the  median  is  limited  at  the  lower  end 
by  the  intersection  of  the  median  with  the  inscribed  circle  of  the  current  front  edge,  (since  this 
is  the  smallest  possible  circumcircle  which  can  be  formed  with  a  triangle  consisting  of  this  edge 
and  a  third  point),  and  at  the  other  extreme  by  the  location  of  the  circumcenter  of  the  current 
waiting  triangle  formed  with  this  edge.  This  extreme  corresponds  to  a  circumcenter  point-insertion 
strategy,  as  described  in  section  3.3,  which  ensures  the  new  point  will  not  fall  coincidentally  close 
to  an  existing  mesh  point.  Figure  19  illustrates  the  mesh  generation  process  at  an  intermediate 
stage  for  the  same  configuration  shown  in  the  previous  section.  The  front  is  clearly  visible,  and  the 
identity  of  the  accepted  and  waiting  triangles  is  fairly  evident.  This  algorithm  uses  the  same  point- 
placement  strategy  as  the  algorithm  of  the  previous  section,  and  thus  should  yield  the  identical  mesh 
as  the  advancing-front  Delaunay  procedure  of  the  previous  section,  provided  the  same  background 
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spacing  function  is  specified  (not  the  case  in  this  example).  The  main  difference  between  these  two 
algorithms  is  that,  in  the  former  case,  the  waiting  triangles  which  are  most  often  highly  skewed,  are 
never  formed,  thus  leaving  an  ungridded  gap  region  between  the  fronts.  This  results  in  the  need  to 
employ  quadtree  data-structures  to  perform  the  search  for  nearby  front  points  and  circumcircles. 
In  the  present  method,  the  domain  is  always  entirely  covered  by  a  triangulation,  which  is  utilized 
in  an  efficient  manner  for  supporting  all  spatial  searching  operations  (implemented  as  neighbor 
walks). 


Figure  19:  Partially  completed  advancing-front  Delaunay  point-insertion  mesh. 

A  notable  feature  of  the  present  methodology  is  that  the  boundary  recovery  procedure  must 
be  performed  immediately  after  the  insertion  of  the  boundary  points,  and  cannot  be  performed  at 
the  end  of  the  mesh  generation  operation,  after  all  interior  points  have  been  added,  as  discussed 
in  section  3.3.  This  is  due  to  the  fact  that  the  boundary  discretization  constitutes  the  initial 
front,  the  location  of  which  is  required  to  guide  the  placement  of  new  points.  Once  the  boundary 
discretization  has  been  recovered,  it  must  be  maintained  throughout  the  insertion  of  all  interior 
mesh  points.  While  in  two  dimensions  this  is  easily  achieved  by  modifying  the  point-insertion 
algorithms  according  to  the  principles  of  constrained  Delaunay  triangulations,  the  lack  of  any 
such  notion  in  three  dimensions  makes  the  situation  more  difficult.  Indeed,  the  initial  mesh  with 
recovered  boundaries  in  general  will  not  represent  a  Delaunay  triangulation,  and  may  lead  to  failure 
of  the  Delaunay-based  point-insertion  algorithms. 

An  alternative  technique  [25,  46]  circumvents  this  issue.  The  approach  is  based  on  relaxing  the 
constraint  that  the  final  mesh  be  as  close  as  possible  to  Delaunay,  and  relies  on  variants  of  the  Green- 
Sibson  point-insertion  algorithm.  The  basic  mesh  generation  procedure  is  similar  to  that  described 
above,  except  that  in  step  8,  a  truncated  Green- Sibson  algorithm  is  employed  to  insert  new^  field 
points.  The  new  point  is  inserted  by  joining  it  to  the  four  vertices  of  its  enclosing  tetrahedron 
(or  five  vertices  of  two  neighboring  tetrahedra  in  the  event  the  point  coincides  with  a  face,  with 
a  similar  extension  for  a  point  coinciding  with  an  edge),  and  then  swapping  faces  and  edges  of 
the  newly  formed  and  subsequently  modified  elements  according  to  the  in-circle  criterion,  until  a 
Delaunay  mesh  is  recovered,  or  the  process  can  no  longer  “improve”  the  mesh.  Boundary  integrity 
is  easily  enforced  in  this  manner,  since  face-edge  swaps  which  alter  the  boundary  discretization 
are  simply  prohibited.  While  the  Bowyer/Watson  algorithm  relies  on  the  Delaunay  property  of 
the  mesh  in  order  to  guarantee  the  convexity  of  the  reconstructed  region  and  the  validity  of  the 
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new  triangulation,  the  Green- Sibson  algorithm  can  still  be  employed  with  initial  constructions 
other  than  Delaunay  tri angulations.  In  such  cases  it  reduces  to  a  local  optimization  technique 
which  attempts  to  approximate  the  Delaunay  triangulation.  By  modifying  the  face-edge  swapping 
criterion,  the  algorithm  can  be  employed  to  construct  other  types  of  triangulations,  such  as  min- 
max  tri  angulations.  (However,  forward  as  well  as  backward  propagation  of  the  swapping  tests 
must  be  implemented  in  such  cases  as  described  in  section  2.1.2).  In  fact,  Delaunay  triangulations 
may  not  be  the  optimal  construction  for  three-dimensional  mesh  generation  problems,  in  spite  of 
aU  the  elegant  properties  associated  with  such  triangulations.  Since  it  is  possible  to  form  a  sliver 
tetrahedron  with  vanishingly  small  volume  for  which  the  circumsphere  is  not  excessively  large  (z.e., 
imagine  four  coplanar  points  on  the  sphere),  Delaunay  triangulations  in  three  dimensions  often 
admit  such  ill-behaved  elements.  This  suggests  that  min-max  triangulations  may  be  better  suited 
for  three-dimensional  grid  generation  purposes. 

On  the  other  hand,  the  use  of  Green-Sibson  type  point-insertion  local-reconnection  algorithms 
for  generating  min-max  triangulations  suffers  from  the  inability  to  recover  the  globally  optimum 
min-max  triangulation.  In  fact,  the  algorithm  typically  gets  caught  very  quickly  in  a  local  optimum 
and  results  in  poor  quality  meshes.  A  key  feature  which  enables  a  closer  approximation  of  the 
global  optimum  and  yields  higher  quality  meshes,  is  the  use  of  an  intermediate  Delaunay  in-circle 
swapping  test,  as  demonstrated  in  [43,  46].  In  this  approach,  each  time  a  new  point  is  inserted,  the 
surrounding  faces  and  edges  are  first  swapped  according  to  the  traditional  Delaunay  in- circle  test, 
and  then  re-swapped  according  to  the  min-max  criterion.  The  use  of  this  intermediate  pseudo- 
Delaunay  construction  serves  to  broaden  the  range  of  influence  of  the  new  point,  thus  providing  a 
more  global  effect,  and  avoiding  the  local  optima  which  typically  plague  the  straight  application  of 
the  min-max  criterion. 

4  Adaptive  Meshing 

Aside  from  the  treatment  of  complex  geometries,  the  second  main  advantage  of  unstructured  meshes 
is  the  ease  with  which  solution-adaptive  meshing  may  be  implemented.  Since  no  inherent  structure 
is  assumed  in  the  representation  of  unstructured  meshes,  mesh  refinement  and  coarsening  may  be 
performed  arbitrarily  in  any  region  of  the  mesh.  The  principle  idea  of  adaptive  meshing  is,  of 
course,  to  enable  a  higher  accuracy  solution  at  lower  cost,  through  a  more  optimal  distribution  of 
grid  points  for  each  computed  solution.  The  entire  procedure  is  akin  to  a  control  problem,  where 
the  output  (the  solution)  drives  the  machinery  which  generates  the  solution  itself.  The  basic  steps 
in  an  adaptive  meshing  solution  strategy  are: 

•  Computation  of  initial  solution. 

•  Estimation  of  local  error  in  solution. 

•  Modification  of  mesh  according  to  estimated  error  values. 

•  Initiahzation  of  solution  on  adapted  mesh. 

•  Resumption  of  numerical  solution  procedure. 

A  complete  adaptive  solution  package  must  therefore  include  a  flow  solution  module,  an  error 
estimation  technique,  and  a  grid  adaptation  (refinement- derefinement)  module,  as  well  as  an  ac¬ 
curate  representation  of  the  actual  configuration  geometry.  Thus,  the  construction  of  an  adaptive 
solution  methodology  represents  alarge  investment  in  software  development,  particularly  in  three 
dimensions. 
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The  need  to  store  geometry  data  is  dictated  by  the  requirements  of  plaeing  new  adaptively 
introduced  boundary  points  on  the  curves  or  surfaces  which  define  the  original  geometry,  rather 
than  simply  on  the  coarse  grid  approximation  to  this  geometry,  as  shown  in  Figure  20. 


Figure  20:  Illustration  of  fine  and  coarse  grid  discretizations  of  curved  boundary. 


The  extraction  and  reduction  of  this  data  from  its  defining  source  (usually  a  CAD  data-package) 
may  prove  to  be  a  difficult  task,  particularly  for  complex  three-dimensional  configurations.  Another 
possibility  is  to  encode  the  geometry  as  a  very  fine  set  of  points  (finer  than  the  highest  expected 
mesh  resolution).  This  seemingly  primitive  technique  may  be  making  a  comeback,  particularly 
with  the  availability  of  rapid  high-resolution  geometry  scanning  devices. 

The  actual  mechanics  for  modifying  the  mesh,  i.e,,  adding,  removing,  and  displacing  mesh 
points,  as  well  as  reconnecting  them,  is  usually  the  topic  which  receives  the  most  elaboration  in 
adaptive  meshing  papers.  It  is  also  perhaps  the  best  understood  of  all  the  required  elements. 
Most  adaptive  meshing  techniques  can  be  implemented  as  extensions  or  enhancements  to  known 
mesh  generation  strategies.  The  formulation  of  an  effective  mesh  refinement  criterion,  on  the  other 
hand,  is  a  difficult  task  which  has  not  been  adequately  resolved  to  date.  The  main  problem  is 
that  an  exact  characterization  of  the  error  requires  a  knowledge  of  the  solution  itself,  which  is 
obviously  impractical.  While  finite-element  error  estimates  have  been  developed  for  simple  effiptic 
problems,  the  difficulty  is  compounded  for  fluid  dynamics  problems  by  the  fact  that  the  governing 
equations  represent  a  coupled  system  of  non-linear  hyperbolic  partial-differential  equations.  A  true 
characterization  of  the  error  would  necessarily  require  information  from  all  flow  variables.  Estimates 
which  rely  on  the  smoothness  properties  of  the  solution  break  down  in  regions  near  discontinuities 
such  as  at  shocks.  Furthermore,  all  error  estimates  rely  on  the  fact  that  the  computed  solution 
is  asymptotically  close  to  the  exact  solution.  The  non-linear  and  hyperbolic  character  of  the 
governing  equations  may,  however,  result  in  situations  wdiere  this  assumption  does  not  hold,  at 
least  locally.  Consider  for  example  a  separation  bubble  which  only  appears  when  certain  upstream 
flow  features  have  been  adequately  resolved,  or  the  diffusion  of  a  wake  profile  due  to  inadequate 
upstream  resolution.  In  such  cases,  the  downstream  solution  is  in  no  w^ay  asymptotically  close  to 
the  exact  solution  in  these  regions. 

The  inadequacy  of  current  error  estimators,  coupled  with  the  fact  that  the  computed  solution 
may  be  far  removed  from  the  exact  solution,  are  the  main  reasons  wdiy  adaptive  meshing  cannot  be 
utilized  in  lieu  of  mesh  generation  itself.  In  other  words,  the  idea  of  initiating  the  calculation  with 
an  extremely  coarse  grid  (possibly  just  fine  enough  to  resolve  the  essential  geometry  topology)  and 
relying  on  adaptive  meshing  to  generate  the  final  mesh,  is  impractical  since  many  flow^  features  will 
never  be  captured  on  such  a  mesh.  The  generation  of  a  suitable  initial  mesh,  with  good  resolution 
in  regions  of  expected  solution  activity,  is  essential  for  good  overall  adaptive  solution  performance. 
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4,1  Refinement  Criteria 


The  most  popular  refinement  criteria  for  fluid-flow  problems  are  essentially  heuristically  derived 
gradient-based  criteria,  which  involve  a  single  or  multiple  physical  flow  variables.  As  an  example, 
the  gradient  of  pressure  can  be  used  to  identify  inviscid  flow  features,  while  the  gradient  of  velocity 
can  be  used  to  track  shear  layers.  In  actual  fact,  it  is  the  undivided  gradient  which  is  employed: 


or,  discretely: 


(2) 


e  =  Au  (.3) 

The  use  of  an  undivided  gradient  ensures  that  the  value  of  e,  which  should  approximate  the  error, 
decreases  as  the  mesh  size  is  reduced.  Since  second-order  methods,  which  are  typically  employed, 
can  represent  hnear  solutions  exactly,  an  undivided  second  difference  may  be  expected  to  yield  a 
better  refinement  criteria: 


c  = 


d'^u 


(4) 


Both  first  and  second  difference  based  refinement  criteria  have  been  tested  and  employed  in  the 
literature.  However,  the  use  of  first  differences  is  more  conservative  (i.e.,  produces  more  refinement  ) 
and  is  often  found  to  work  better  in  practice.  Lohner  [47]  advocates  the  use  of  a  non-dimensional 
criterion  which  is  designed  to  assign  equal  weight  to  weak  as  well  as  strong  flow  features.  This 
is  achieved  by  forming  the  ratio  of  second  and  first  derivatives.  The  refinement  indicator  can  be 
written  as: 


€ 


I  I 
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where  u  represents  some  neighborhood  average  of  the  flow  variable  u,  and  a  is  a  small  parameter. 
This  extra  term  in  the  denominator  acts  as  a  “noise”  filter  which  avoids  triggering  refinement  in 
regions  of  small  solution  oscillations.  This  type  of  refinement  criterion  has  been  used  extensively 
for  two-  and  three-dimensional  transient  flow  solutions  involving  shock  waves. 

An  interesting  study  of  the  effectiveness  of  various  refinement  criteria  for  steady-state  problems 
can  be  found  in  the  paper  by  Warren  et  al.  [48].  The  error  levels  in  various  adapted  mesh  solutions 
using  different  refinement  criteria  were  assessed  by  comparing  the  adapted  mesh  solutions  to  the 
“exact  ’  solution  (i.e.,  a  solution  computed  on  a  highly  resolved,  globally  refined  mesh).  Their 
experiments  indicate  that  the  modification  of  the  undivided  differences  in  equation  (3)  to  include 
a  local  mesh  length-scale  such  as: 


€  =  Au  •  A.t  (6) 

produces  a  more  effective  refinement  criterion.  This  is  partially  due  to  the  fact  that  while  the  simple 
undivided  difference  form  of  equation  (3)  decreases  in  magnitude  as  the  mesh  is  refined  in  smooth 
regions  of  flow,  it  remains  approximately  constant  in  the  vicinity  of  shock  waves,  since  the  shock 
wave  profile  steepens  as  the  mesh  is  refined,  and  the  jumps  remain  almost  constant.  However,  even 
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in  regions  of  smooth  flow,  the  additional  length  scale  weights  large  cells  more  heavily  than  small 
cells,  and  drives  the  adaptation  process  closer  towards  global  refinement.  The  introduction  of  an 
extra  length  scale  thus  produces  a  more  conservative  criterion  (f.e.,  results  in  more  refinement) 
and  thus  greater  reduction  of  the  error  can  be  expected.  This  asymptotic  tendency  towards  global 
refinement  is  actually  the  behavior  which  is  desired  for  steady-state  problems.  While  adaptive 
meshing  is  typically  thought  of  as  a  process  that  refines  only  local  regions  of  the  mesh,  this  is 
characteristic  of  the  initial  stages  of  an  adaptive  mesh  procedure,  where  the  appropriate  distribution 
of  mesh  resolution  is  set  up  to  match  the  solution.  Once  this  has  been  achieved,  and  the  error  is 
presumably  equipartitioned  throughout  the  domain,  any  increase  in  solution  accuracy  can  only 
come  about  from  global  refinement.  Refinement  criteria  must  therefore  be  able  to  reproduce  this 
behavior. 

Although  gradient-based  refinement  criteria  have  been  employed  successfully  in  computational 
fluid  dynamics,  such  techniques  are  not  well  founded  and  are  far  from  optimal.  The  effectiveness 
of  a  refinement  criterion  is  measured  in  terms  of  accuracy  delivered  for  the  number  of  mesh  points 
required.  In  the  above  discussion,  the  more  conservative  gradient-based  criteria  are  obviously  the 
most  successful  at  reducing  the  solution  error,  but  also  result  in  excessive  refinement,  and  are  not 
necessarily  the  most  efficient.  The  construction  of  more  optimal  refinement  criteria  will  require 
the  development  of  better  error  estimation  techniques  most  likely  through  the  use  of  some  type  of 
extrapolation  method  [49]. 

4.2  Mesh  Adaptation  Techniques 

Mesh  adaptation  may  involve  the  addition  of  extra,  vertices,  the  removal  of  vertices,  the  redistri¬ 
bution  of  existing  vertices  and  the  reconnection  of  mesh  vertices.  The  character  of  the  problem 
to  be  solved  dictates  the  requirements  of  the  mesh  adaptation  strategy.  For  example,  steady-state 
problems  usually  involve  a  small  number  of  adaptation  phases  as  part  of  a  lengthy  solution  process. 
Therefore,  relatively  sophisticated  (i.e..  more  optimal)  adaptation  strategies  can  be  employed,  such 
as,  in  the  extreme  case,  complete  mesh  regeneration.  Mesh  refinement  procedures  are  most  impor¬ 
tant  here,  while  de-refinement  has  only  a  minor  effect  and  can  often  be  omitted  for  steady-state 
cases.  For  transient  problems,  mesh  adaptation  must  be  performed  every  several  time-steps,  and 
thus  efficiency  is  much  more  important  than  optimality.  Mesh  refinement  and  de-refinement  are 
both  essential  for  transient  cases.  Furthermore,  the  accuracy  of  interpolation  from  the  original 
mesh  to  the  refined  mesh  affects  the  solution  accuracy  (unlike  the  steady-state  case),  and  thus 
accurate  transfer  schemes  are  required.  These  requirements  have  often  lead  to  the  use  of  simple 
element  subdivision  schemes  for  transient  flows. 

For  steady-state  as  well  as  transient  problems,  mesh  adaptation  strategies  are  generally  based 
on  one  of  the  previously  discussed  mesh  generation  procedures,  and  can  be  implemented  as  an 
extension  of  the  originating  procedure. 

4.2.1  Delaunay  Point-Insertion  Based  on  Solution  Gradients 

The  Bowyer/Watson  and  Green-Sibson  algorithms  described  in  section  2.1  enable  the  introduction 
of  new  points  in  any  regions  of  an  existing  mesh,  and  are  thus  natural  candidates  for  mesh  adap¬ 
tation  algorithms.  A  particular  implementation  of  an  adaptive  mesh  solution  strategy  using  the 
Bowyer/Watson  algorithm  performed  by  the  author  [-50]  is  described  below.  Assuming  an  initial 
Delaunay  triangulation  mesh  has  been  constructed,  the  flow  solution  is  computed  on  this  mesh. 
Using  the  undivided  gradient  of  density  as  a  refinement  criterion  {c.f.  equation  (3)),  each  edge  of 
the  mesh  is  then  examined  and  flagged  for  refinement  if  the  difference  of  density  along  the  edge 
is  larger  than  some  threshold  value,  which  is  set  proportional  to  the  average  of  all  density  dif¬ 
ferences  taken  over  all  mesh  edges.  This  simple  edge-based  refinement  strategy  tends  to  produce 
non-isotropic  point  distributions  which  triangulate  poorly.  Therefore,  a  more  isotropic  refinement 
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strategy  is  desired.  This  is  achieved  by  effecting  a  single  loop  over  the  mesh  triangles  and  flagging 
the  remaining  edges  for  refinement  of  any  triangle  which  has  a  single  or  two  refined  edges.  Note  that 
this  operation  requires  the  use  of  a  cell-to-edge  data-structure,  in  addition  to  the  ceU-to-node  data- 
structure  employed  by  the  Delaunay  triangulation  algorithm,  and  the  edge-to-node  data-structure 
employed  by  the  flow  solution  algorithm. 

For  each  flagged  edge,  a  new  point  is  created  at  the  midpoint  of  the  edge.  For  boundary  edges, 
these  points  are  repositioned  onto  the  spMne  curve  which  defines  the  original  geometry,  as  depicted 
in  Figure  20.  The  new  points  are  put  in  a  list,  and  then  inserted  and  triangulated  sequentially  into 
the  mesh  using  the  Bowyer/Watson  algorithm.  The  new  refined  mesh  is  smoothed  by  repositioning 
the  points  according  to  a  Laplacian  filtering  technique.  The  solution  is  then  interpolated  in  a 
piecewise  linear  manner  from  the  original  mesh  to  the  new  finer  mesh  (using  the  interpolation 
operators  constructed  for  the  multigrid  algorithm  as  described  in  the  next  chapter)  and  the  flow  is 
solved  on  the  new  mesh.  The  entire  procedure  is  repeated  several  times,  until  the  desired  accuracy 
(or  grid  size)  is  attained. 


Figure  21:  Adaptively  generated 
mesh  for  NACA  0012  airfoil. 


Figure  22:  Adaptively  generated 
mesh  for  idealized  four  element  airfoil 
geometry. 


Figure  21  illustrates  the  final  mesh  obtained  by  this  procedure,  after  three  adaptive  refinement 
passes,  for  the  computation  of  transonic  flow  over  a  NACA  0012  airfoil.  The  flow  features  are 
rather  well  defined  in  this  case,  and  the  various  levels  of  refinement  are  evident  in  the  final  mesh. 
Figure  22  illustrates  an  adaptive  mesh  obtained  for  the  computation  of  inviscid  subsonic  flow  over 
a  multi-element  airfoil.  Here,  the  flow  is  smooth  and  features  are  not  weU  defined.  The  various 
levels  used  to  construct  the  adapted  mesh  are  not  as  evident,  and  the  procedure  results  in  a  smooth 
variation  of  the  mesh  in  most  regions  of  the  domain.  This  same  procedure  extends  readily  to  three 
dimensions.  Figure  23  depicts  an  adapted  tetrahedral  mesh  employed  to  compute  the  transonic 
flow  over  an  ONER  A  M6  wing.  This  mesh  contains  a  total  of  173,412  points,  and  was  obtained 
through  two  applications  of  adaptive  refinement.  The  double  shock  structure  in  the  flow  solution 
is  clearly  reflected  in  the  mesh  refinement  pattern. 


Figure  23:  Adaptively  generated  mesh  in  three  dimensions  for  ONERA  M6  wing. 


4.2.2  Extensions  to  Steiner  Triangulation 

Any  mesh  generation  technique  which  relies  on  the  use  of  a  background  element-size  distribution 
function  can  be  extended  in  a  natural  manner  for  adaptive  meshing  problems.  Consider  the  Steiner 
triangulation  methods  described  in  section  3.3,  where  new  points  are  continually  added  to  the  mesh, 
until  elements  which  satisfy  the  background  size-function  are  obtained.  Once  the  flow  solution  has 
been  obtained  on  the  initial  mesh,  the  background  function  may  be  modified  to  reflect  the  result 
of  the  application  of  a  refinement  criterion  to  this  solution,  (be.,  the  background  function  can  be 
modified  to  specify  small  element  sizes  in  regions  of  high  flow  gradients  or  solution  error).  A  new 
adapted  mesh  may  then  be  constructed  by  simply  restarting  the  Steiner  triangulation  algorithm, 
and  resuming  triangulation  until  convergence  to  the  new  background  function  is  achieved.  This 
strategy,  of  course,  only  produces  additional  refinement.  In  principle,  mesh  de- refinement  may 
also  be  achieved  by  removing  points  sequentially,  and  retriangulating  the  resulting  convex  cavity 
for  each  deleted  point,  in  regions  where  the  mesh  element-size  is  smaller  than  that  prescribed  by 
the  background  function.  This  may  be  achieved,  for  example,  by  applying  the  Tanemura-Merriam 
Delaunay  triangulation  algorithm  to  the  vertices  of  the  cavity  which  is  created  upon  the  removal  of  a 
point  from  the  current  mesh.  In  three  dimensions,  however,  it  is  possible  to  encounter  cavities  which 
cannot  be  tetrahedralized,  and  a  more  sophisticated  approach  may  be  required.  An  alternative  is 
to  maintain  a  data-structure  which  encodes  a  hierarchical  history  of  the  triangulation  construction, 
as  suggested  by  Barth  [51]. 

4.2.3  Adaptive  Remeshing 

The  background  element-size  distribution  function  of  traditional  advancing-front  methods  may 
also  be  modified  in  a  solution  adaptive  manner.  For  background  functions  defined  on  a  (coarse) 
unstructured  grid,  one  approach  is  to  define  the  new  function  on  the  current  mesh  slated  for 
refinement  [26].  The  advancing-front  mesh  generation  software  can  then  be  utilized  to  regenerate 
a  new  mesh  from  scratch,  using  the  newly  defined  background  function  to  determine  the  solution- 
adaptive  mesh  distribution.  A  less  expensive  alternative  is  to  use  local  remeshing  [34].  In  this 
approach,  the  current  mesh  is  removed  in  various  regions  of  the  domain  where  the  discrepancy 
between  current  and  desired  mesh  resolution  is  large.  New  local  meshes  are  then  generated  in 
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these  void  regions  using  the  advancing-front  method,  where  the  initial  front  is  determined  by  the 
boundary  between  meshed  and  void  regions. 

4.2.4  h-Refinement  or  Subdivision  Techniques 

h-refinement  techniques,  which  rely  on  the  subdivision  of  mesh  elements  according  to  a  prede¬ 
termined  set  of  primitive  operations,  are  both  simple  and  efficient.  These  methods  result  in  fine 
grid  elements  which  are  fully  nested  with  their  forming  coarse  grid  elements,  as  well  as  fine  grids 
which  contain,  as  a  subset  of  their  vertices,  all  the  coarse  grid  vertices.  These  properties  enable 
a  very  accurate  and  efficient  transfer  of  variables  from  one  grid  to  another,  as  well  as  a  simple 
framework  for  encoding  history  effects  to  determine  parent-child  relationships  between  cells.  This, 
in  turn,  enables  efficient  use  of  de-refinement  techniques.  Efficiency,  de- refinement,  and  accuracy 
of  interpolation  make  these  methods  ideal  for  transient  problems.  The  fully  nested  property  of 
these  techniques  has  also  been  exploited  for  constructing  unstructured  multigrid  algorithms  based 
on  such  meshes  (see  chapter  on  multigrid  methods). 

The  main  drawback  of  such  methods  is  the  possibility  of  generating  ill-shaped  elements  and 
meshes  with  arbitrarily  high  connectivity  at  isolated  vertices.  In  order  to  avoid  such  situations, 
strict  rules  on  the  permitted  element  subdivision  types  must  be  enforced.  A  set  of  primitives  for 
the  subdivision  and  de-refinement  of  three-dimensional  tetrahedra  are  given  by  Lohner  [47],  and 
reproduced  in  Figure  24. 


Figure  24:  Rules  for  subdivision  refinement  of  three  dimensional  tetrahedral  meshes. 


1:2,  1:4,  and  1:8  refinements  are  permitted.  Subsequent  refinements  of  2:4,  2:8,  and  4:8  are  also 
permitted.  However,  if  a  1:2,  1:4,  or  2:4  configuration  is  to  be  refined  further,  it  must  first  be 
transformed  to  a  1:8  refinement,  and  then  its  child  elements  refined  recursively.  When  applying 
these  rules  to  a  set  of  elements  in  a  mesh,  care  must  be  taken  to  ensure  compatible  refinement 
patterns  are  obtained  on  neighboring  elements,  and  to  avoid  jumps  larger  than  8:1  between  neigh¬ 
boring  elements.  This  is  achieved  by  iteratively  modifying  the  refinement  types  on  groups  of  mesh 
elements  until  a  compatible  pattern  is  achieved. 
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Figure  25:  Example  of  three-dimensional  mesh  adaptation  by  subdivision  for  transient 
shock-wave  problem  over  an  unfriendly  configuration.  (Reproduced  from  [52]  with  permis¬ 
sion). 

Figure  25  illustrates  an  example  of  a  transient  shock  wave  problem  computed  with  an  adaptive 
h-refinement  technique,  taken  from  [52],  using  the  technique  described  above.  Another  method  for 
avoiding  ill-shaped  elements  and  poor  mesh  connectivity  is  to  perform  an  edge-swapping  or  face- 
edge  swapping  operation  in  two  dimensions  or  three  dimensions,  respectively,  after  the  element 
subdivision  operations,  in  order  to  improve  the  mesh  topology.  While  substantial  improvement 
in  mesh  quality  can  be  obtained,  particularly  if  the  strict  refinement  rules  have  been  relaxed, 
the  nested  property  of  the  meshes  is  lost,  and  the  parent-child  history  essential  for  de-refinement 
becomes  more  complicated  to  recover. 

4.2,5  Mesh  Movement  Techniques 

Mesh  movement  techniques  involve  the  movement  of  mesh  vertices  in  order  to  improve  the  dis¬ 
cretization  accuracy,  by  drawing  points  towards  regions  of  high  solution  error,  and  away  from 
regions  of  excessive  resolution.  These  are  often  performed  without  altering  the  mesh  connectivity, 
which  makes  then  particularly  simple  to  implement.  While  they  have  proved  useful  for  structured 
meshes  (where  the  mesh  structure  must  be  preserved),  they  are  rarely  used  in  this  form  for  un¬ 
structured  meshes,  since  mesh  refinement  is  easily  implemented.  They  are  most  often  employed  in 
conjunction  with  other  adaptive  meshing  techniques  as  a  smoother,  in  order  to  improve  local  grid 
quality.  A  simple  construction  of  Laplacian-type  smoothing  operations  for  two-dimensional  meshes 
is  given  by  updating  the  coordinates  of  each  mesh  point  as: 


where  the  summation  is  over  all  the  neighboring  vertices  of  the  considered  point.  This  procedure 
mimics  a  Jacobi  iteration  for  a  Laplacian  operator  on  the  coordinates  of  the  grid  points,  and  can  be 
applied  iteratively  for  all  grid  points.  This  particular  formulation  does  not  exclude  the  possibility 
of  forming  negative  area  elements.  Although  formulations  have  been  developed  for  excluding  such 
possibilities,  they  are  often  expensive  and  difficult  to  solve.  A  simpler  approach  is  to  unsmooth 
grid  points  in  regions  where  negative  elements  are  created.  This  approach,  combined  with  the  use 
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of  face-edge  swapping  for  improving  grid  connectivity,  can  substantially  improve  mesh  quality  (c./. 
Figure  14). 

While  mesh  movement  methods  have  seldom  been  employed  as  a  solution-adaptive  meshing 
procedure  for  unstructured  meshes,  a  notable  exception  is  the  use  of  these  techniques  for  shock 
fitting  purposes  [53,  54].  If  the  mesh  points  of  a  triangular  mesh  are  displaced  in  such  a  manner  as 
to  align  the  mesh  edges  with  the  shocks  in  the  flow,  then  it  is  possible  to  capture  the  shock  over  a 
single  mesh  cell  interface  using,  for  example,  a  cell-centered  scheme  with  a  Rieman  interface  solver. 


Figure  26:  Original  mesh,  adapted  mesh,  and  corresponding  solutions,  using  mesh- 
movement  shock  fitting  approach.  (Reproduced  from  [54]  with  permission). 


Figure  26  illustrates  the  use  of  this  technique  for  capturing  a  bow  and  fishtail  shock  about  a 
supersonic  airfoil.  Not  only  can  very  weak  shocks  be  weU  captured,  but  the  required  number  of 
grid  points  and  computational  effort  is  much  lower  than  for  mesh  refinement  techniques. 

5  Stretched-Mesh  Generation 

The  drive  towards  full  Navier- Stokes  solvers  has  necessitated  the  development  of  stretched  grid 
generation  techniques  in  order  to  resolve  the  thin  boundary-layer  and  wake  regions  which  are 
characteristic  of  high-Reynolds-number  viscous  flows.  Proper  boundary-layer  and  wake  resolution 
usually  requires  mesh  spacing  several  orders  of  magnitude  smaller  in  the  direction  normal  to  the 
boundaries  than  in  the  streamwise  direction,  resulting  in  large  cell  aspect-ratios  in  these  regions. 
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The  application  of  the  triangulation  methods  described  in  the  preceding  sections  to  these  types 
of  mesh-point  distributions  results  in  awkward  triangulations,  which  are  ill-suited  for  numerical 
computations.  This  is  due  to  the  fact  that  most  of  these  methods  have  been  conceived  for  the 
generation  of  isotropic  meshes,  and  are  ill-suited  for  the  generation  of  high-aspect-ratio  elements. 
Consider,  for  example,  the  Delaunay  triangulation  of  such  a  point-set.  In  two  dimensions,  Delaunay 
triangulations  are  equivalent  to  max-min  triangulations,  z.e.,  they  correspond  to  the  triangulation 
which  maximizes  the  smallest  angles  of  all  elements.  Since  high  aspect-ratio  triangles  necessar¬ 
ily  contain  at  least  one  very  small  angle,  Delaunay  triangulations  actually  attempt  to  avoid  the 
generation  of  such  elements. 

There  are  two  possible  approaches  in  dealing  with  these  difficulties.  The  first  is  to  make  use  of 
alternative  element  types  in  regions  of  high  mesh  stretching,  such  as  quadrilaterals  in  two  dimensions 
and  prisms  or  hexahedra  in  three  dimensions  [55,  56].  The  second  approach  is  to  define  the  types 
of  triangle/tetrahedral  elements  which  are  desirable  for  stretched  mesh  generation,  and  to  modify 
existing  methods  or  devise  new  techniques  for  generating  meshes  which  contain  such  elements. 

In  practice,  both  approaches  are  viable.  In  fact,  the  two  approaches  are  less  distinct  than  may 
be  expected,  since  it  is  always  possible  to  subdivide  a  mixed-element  mesh  into  a  fully  triangular 
or  tetrahedral  mesh,  and  similarly  most  highly  stretched  triangular/tetrahedral  meshes  can  be 
transformed  into  mixed  quadrilateral/triangular  or  prismatic-tetrahedral  meshes  by  identifying 
appropriate  mesh  edges  for  removal  [15,  43].  In  fact,  the  definition  of  a  stretching  direction  in 
itself  implies  a  certain  degree  of  local  structure  in  the  mesh,  which  enables  the  simple  extraction  of 
quadrilaterals/prisms  from  a  triangular /tetrahedral  mesh.  Similarly,  no  loss  of  flexibility  is  implied 
by  the  use  of  semi-structured  quadrilateral  or  prismatic  meshes  in  highly-stretched  regions.  The 
main  drawbacks  of  using  mixed-element  meshes  is  the  lack  of  homogeneity  in  the  grid  structure, 
which  may  complicate  procedures  such  as  flow  solution  and  mesh  adaptation.  The  advantages  are 
reduced  overheads  for  these  same  procedures,  due  to  the  lower  overall  connectivity  of  these  mesh 
elements  and  the  resulting  reduction  in  the  number  of  mesh  edges.  In  three  dimensions,  the  use 
of  prismatic  elements  are  particularly  attractive,  since  this  permits  the  use  of  triangular  surface 
meshes. 

If  fully  triangular /tetrahedral  meshes  are  to  be  employed,  a  characterization  of  the  optimal 
stretched  triangle  shape  needs  to  be  defined.  In  [57],  it  is  showui  how  the  accuracy  of  a  two- 
dimensional  finite-element  approximation  on  triangular  elements  degrades  as  the  maximum  angle 
of  the  elements  increases. 


(a) 


(b) 


Figure  27:  Illustration  of  two  types  of  high  aspect  ratio  triangles:  non-obtuse  (a)  and 
obtuse  (b). 


Thus,  as  shown  in  Figure  27,  stretched  obtuse  triangles  which  contain  one  large  angle  and  two  small 
angles  are  to  be  avoided,  while  stretched  nearly  right-angle  triangles,  with  one  small  angle  and  two 
nearly  right  angles  are  to  be  preferred,  (Similarly,  one  may  infer  that  the  types  of  tetrahedra  which 
result  from  the  subdivision  of  a  thin  prism  in  3D  are  desirable).  Thus,  a  triangulation  procedure 
which  avoids  obtuse  triangles  in  favor  of  right-angle  triangles  would  be  desirable.  It  should  be 
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noted  ,  however,  that  even  if  a  “perfect”  triangulation  scheme  exists,  acceptable  mesh  elements  will 
only  be  obtained  provided  the  vertices  are  positioned  appropriately.  Thus,  most  stretched-mesh 
generation  techniques  involve  careful  point-placement  strategies  as  well  as  modified  triangulation 
schemes. 

5.1  Stretched  Delaunay  Point-Insertion 

One  of  the  earliest  techniques  for  generating  highly  stretched  triangular  meshes  [36,  37,  58]  makes 
use  of  a  locally-mapped  Delaunay  triangulation  point-insertion  method.  The  basic  idea  is  to  de¬ 
fine  a  local  stretching  vector  (direction  and  magnitude  of  stretching)  at  each  vertex.  New  points 
are  inserted  into  the  triangulation  using  the  Bowyer/Watson  Delaunay  triangulation  algorithm. 
However,  this  triangulation  is  carried  out  in  a  transformed  space,  which  is  obtained  by  locally 
remapping  the  physical  space  according  to  the  values  of  the  stretching  vector  in  the  vicinity  of 
the  new  point.  This  results  in  a  nearly  isotropic  Delaunay  triangulation  in  the  transformed  space, 
but  in  a  stretched  Delaunay  triangulation  in  physical  space,  where  triangle  circumcircles  become 
circumeUipses,  as  shown  in  Figure  28. 


Figure  28:  CircumeUipse  of  a  stretched  Delaunay  triangulation. 


The  problem  of  constructing  a  global  mapping  for  complex  domains  with  variable  distribntions  of 
stretching  represents  a  formidable  task.  However,  since  the  the  Delaunay  triangulation  is  a  local 
construction,  the  problem  can  be  formulated  in  a  simpler  manner,  using  local  mappings.  Each 
time  a  new  point  is  to  be  inserted,  a  local  mapping  may  be  constructed,  using  the  local  values 
of  the  stretching  vectors,  which  may  be  averaged  and  taken  as  constant  throughout  the  region 
to  be  restructured.  This  enables  the  use  of  a  variable  stretching  distribution,  but  requires  the 
distribution  of  stretching  to  be  smooth,  and  to  vary  slowly  with  respect  to  the  point-distribution. 
In  fact,  a  close  coupling  between  the  stretching  and  mesh-point  distributions  is  required  in  order  to 
ensure  the  generation  of  appropriately  shaped  triangular  elements.  In  the  present  approach,  this 
is  achieved  by  generating  both  mesh-point  and  stretching  distributions  from  a  set  of  overlapping 
stretched  structured  meshes.  The  original  method  [37]  can  be  summarized  as  follows: 


Step  1:  Generate  a  highly  stretched  structured  mesh  about  each  geometry  component. 

Step  2:  Filter  out  far-held  and  downstream  wake  points  in  these  structured  meshes. 

Step  3:  Dehne  a  stretching  vector  at  each  remaining  point,  based  on  the  local  structured  grid-ceU 
orientation,  stretching  direction,  and  aspect-ratio  (stretching  magnitude). 

Step  4:  Using  the  Bowyer/Watson  algorithm,  construct  the  (unstretched)  Delaunay  triangulation 
of  this  set  of  points. 

Step  5:  Smooth  the  distribution  of  stretching  vectors  by  performing  several  passes  of  averaging 
stretchings  with  their  neighboring  values. 
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Step  6:  Swap  the  edges  of  the  mesh  according  to  the  Delaunay  in  circle  criterion,  measured  in  the 
locally  mapped  space,  as  defined  by  the  average  local  stretching  vector  (i.e.,  a  stretched  Delaunay 
in  ellipse  criterion). 

Step  7:  Smooth  the  mesh-point  distribution  and  reswap  the  edges  (this  step  may  be  repeated 
several  times). 

The  use  of  point-sets  derived  from  structured  meshes  ensures  that  the  points  will  be  distributed  in 
a  configuration  which  favors  the  formation  of  nearly  right-angle  triangles  in  the  highly  stretched 
regions  of  the  mesh,,  such  as  near  the  walls  and  in  the  wake  regions.  Step  .3  ensures  the  compatibility 
between  the  stretching  distribution  (and  thus  triangulation  criterion)  and  point  distribution,  while 
step  5  is  necessary  to  guarantee  a  smooth  distribution  of  stretching  in  regions  where  the  structured 
meshes  overlap.  It  is  noteworthy  that  an  initial  triangulation  is  required  in  order  to  accomplish 
this  smoothing.  Figure  29  illustrates  a  stretched  unstructured  mesh  produced  by  this  technique  for 
a  simple  airfoil  configuration. 


Figure  29:  Stretched  unstructured  mesh  produced  by  mapped  Delaunay  trlangulation 
procedure  for  RAE  2822  airfoil. 

The  present  method  permits  a  smooth  transition  between  stretched  and  unstretched  regions,  as 
well  as  between  two  neighboring  stretched  regions,  and  enables  a  straight-forward  implementation 
of  adaptive  meshing  techniques.  New  mesh  points  may  be  inserted  into  the  existing  mesh  by  first 
assigning  them  a  stretching  vector  taken  as  the  average  of  the  neighboring  stretching  vectors,  and 
then  using  the  Bowyer/Watson  algorithm  in  the  locally  stretched  space  to  triangulate  the  new 
point.  The  final  mesh  may  be  post-processed  with  several  passes  of  smoothing  and  edge-swapping. 
An  example  of  an  adaptively  generated  stretched  unstructured  mesh  about  a  four-element  airfoil 
is  illustrated  in  Figure  30. 


36 


Figure  30:  Adapted  stretched  mesh  for  computing  viscous  flow  over  four-element  airfoil. 

In  retrospect,  this  mesh  generation  approach  can  be  viewed  as  the  construction  of  an  initial  (De¬ 
launay)  triangulation,  followed  by  an  edge-swapping  reconnection  phase  governed  by  a  criterion 
determined  by  a  distribution  of  stretching.  In  this  respect,  it  is  similar  to  some  of  the  more  recent 
methods  for  constructing  stretched  triangulations. 

5.2  Hybrid  Methods 

One  of  the  most  evident  techniques  for  constructing  stretched  meshes  for  viscous  flows  is  to  use  a 
hybrid  technique  which  makes  use  of  a  structured  mesh  in  the  thin  boundary-layer  and  wake  regions, 
which  require  high  degrees  of  stretching,  and  an  unstructured  mesh  in  the  regions  of  inviscid  flow, 
where  isotropic  constructions  are  desirable.  In  two  dimensions,  the  structured  portion  of  the  mesh 
consists  of  quadrilaterals,  while  in  three  dimensions  either  hexahedra  or  prismatic  (semi-structured) 
elements  may  be  used.  These  elements  may  be  later  subdivided  into  triangles  or  tetrahedra  to 
produce  a  homogeneous  mesh,  or  alternatively  they  may  be  retained,  and  the  flow  solver  may  be 
modified  to  capitalize  on  the  local  grid  structure. 

An  example  of  a  hybrid  approach  to  stretched  mesh  generation  consists  of  generating  a  local 
structured  mesh  using  a  hyperbolic  mesh  generation  technique,  up  to  a  prespecified  distance  away 
from  the  boundary  surfaces,  and  an  advancing-front  unstructured  mesh  generation  technique  to 
complete  the  mesh  by  filling  in  the  remainder  of  the  domain  with  an  isotropic  unstructured  mesh 
[56].  The  problems  associated  with  hybrid  mesh  constructions  involve  the  adequate  definition  of 
the  thickness  of  the  structured  mesh  layer,  as  well  as  difficulties  involved  in  concave  regions  or 
regions  where  neighboring  boundaries  are  in  close  proximity,  which  may  result  in  overlapping  of 


the  structured  meshes.  Another  difficulty  consists  of  producing  a  mesh  with  a  smooth  transition 
between  stretched  structured  regions  and  isotropic  unstructured  regions.  Many  of  these  difficulties 
can  be  overcome  by  employing  a  structured  grid  layer  of  variable  thickness. 

5.3  Semi-Structured-Unstructured  Hybrid  Meshes 

The  hybrid  structured-unstructured  approach  may  be  made  more  flexible  by  employing  so-called 
semi-structured  meshes  in  the  stretched  mesh  regions  [59].  Semi-structured  meshes  are  essentially 
structured  meshes  with  a  variable  thickness  or  normal  resolution.  A  semi-structured  mesh  may 
be  constructed  by  first  generating  a  structured  mesh  about  each  geometry  component,  over  a 
distance  which  covers  a  significant  portion  of  the  domain  a  distance  much  greater  than  the 

boundary-layer  thickness),  and  then  removing  the  regions  of  the  mesh  which  overlap  with  neighbor¬ 
ing  structured  meshes  and/or  geometry  components,  as  well  as  regions  of  the  mesh  where  stretching 
is  not  required,  z.e.,  where  the  structured  mesh  cell  aspect-ratios  are  less  than  unity.  The  remaining 
portion  of  the  domain  is  then  gridded  with  an  isotropic  triangulation  scheme.  This  results  is  an 
automatic  determination  of  the  inner-mesh  thickness,  and  a  more  smooth  transition  between  the 
inner  and  outer  meshes. 


Figure  31:  Semi-structured  mesh  gener-  Figure  32:  Final  Unstructured  mesh  gen- 

ated  prior  to  completion  of  isotropic  region.  erated  using  semi-structured  approach. 

An  example  of  a  semi-structured  mesh  generated  by  the  author  is  depicted  in  Figures  31  and 
32.  The  outer  boundary  is  somewhat  ragged,  but  this  presents  no  difficulty  for  the  unstructured 
mesh  approach  which  is  used  to  complete  the  mesh.  In  this  case,  the  advanciiig-front  Delaunay 
triangulation  scheme  was  employed  in  the  isotropic  regions  of  the  domain,  and  elements  of  the  semi- 
structured  mesh  have  been  subdivided  into  triangles.  A  noteworthy  point  concerns  the  structure 
of  the  mesh  in  the  wake  regions.  While  highly  stretched  meshes  are  required  in  wake  regions, 
unless  these  are  to  be  carried  to  the  outer  boundary  of  the  domain,  they  must  be  blended  with 
the  isotropic  portion  of  the  mesh  in  a  smooth  manner.  In  the  present  implementation,  this  is 
accomplished  by  coarsening  the  wake-point  distribution  in  the  normal  direction  increasingly  in  the 
downstream  direction,  and  readjusting  the  connectivity.  This  operation  is  performed  prior  to  the 
generation  of  the  isotropic  portion  of  the  mesh,  since  it  results  in  a  modification  of  the  initial 
boundary  discretization  of  the  remaining  area  to  be  gridded  isotropically. 

The  drawbacks  of  this  method  are  the  requirement  of  being  able  to  generate  structured  meshes 
about  arbitrary  geometry  components,  the  task  of  locating  overlapping  structured  mesh  regions, 
and  the  ability  of  the  method  to  merge  two  neighboring  semi-structured  meshes. 
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5.4  Advancing-Layers  Method 

The  advandng-layers  method  [60,  61,  62,  63,  64]  represents  a  generalization  of  the  concept  of  semi- 
structured  meshes.  Rather  than  construct  an  inner  semi-structured  mesh  by  truncating  a  structured 
mesh,  the  advancing-layers  technique  uses  structured  hyperbolic  mesh  generation  principles  locally 
at  each  boundary  edge  or  face  to  advance  a  new  layer  in  the  normal  direction,  thus  creating  a  new 
ceU.  Advancing-layers  can  also  be  viewed  as  a  modification  to  the  unstructured  advancing-front 
algorithm  to  create  stretched  quadrilateral  (or  prismatic)  elements  using  different  point-placement 
and  reconnection  strategies.  This  is  achieved  by  capitalizing  on  the  similarities  between  hyperbolic 
mesh  generation  and  advancing-front  techniques. 

The  point  placement  strategy  of  the  advancing-layers  method  relies  on  the  surface  normals  of 
the  boundary  discretization.  A  line  segment  normal  to  the  boundary  surface  is  first  associated  with 
each  boundary  point.  These  normals  are  often  smoothed  by  averaging  their  orientations  with  those 
of  their  neighbors.  Progressive  smoothing,  which  vanishes  near  the  boundary  and  increases  away 
from  the  boundary  is  typicaUy  employed.  Quadrilateral  elements  (or  prismatic  elements  in  three 
dimensions)  are  created  by  placing  new  points  along  the  boundary  normals  and  connecting  them 
to  neighboring  points  in  a  specified  pattern. 

The  layers  are  advanced  out  from  the  original  discretized  boundaries  with  an  increasing  step- 
size  which  results  in  progressively  decreased  stretching.  For  each  normal  stack  of  cells  associated 
with  a  boundary  edge  or  face,  the  advancing  process  is  terminated  either  when  an  opposing  front 
is  encountered,  or  when  the  cell  aspect-ratio  becomes  close  to  unity.  Once  the  advancing-layers 
procedure  terminates,  the  remainder  of  the  domain  is  filled  in  with  an  isotropic  unstructured  mesh, 
as  in  the  previous  methods.  While  the  rate  of  growth  of  the  advancing  layers  is  typically  prescribed, 
a  background  function  is  used  to  determine  the  resolution  of  the  initial  boundary  discretization 
as  well  as  to  control  that  of  the  unstructured  portion  of  the  mesh.  When  the  advancing-layers 
phase  terminates,  the  existing  mesh  looks  much  like  the  semi-structured  or  truncated  hyperbolic 
structured  mesh  of  Figure  31.  A  three-dimensional  example  of  the  advancing  layers  method,  taken 
from  [64],  is  depicted  in  Figure  33,  for  a  four-element  wing  configuration,  showing  the  surface  mesh, 
and  the  partially  completed  mesh  of  stretched  layers.  As  in  the  semi-structured  mesh  procedure, 
the  layers  may  be  conserved  as  quadrilaterals  or  prisms  in  two  or  three  dimensions  respectively,  or 
they  may  be  subdivided  into  triangles  or  tetrahedra. 


Figure  33:  Illustration  of  advancing-layers  method  for  three-dimensional  unstructured 
mesh  generation  about  segmented  wing  geometry.  (Reproduced  from  [64]  with  permission). 
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5.5  Advancing-Front  Min-Max  Tidangulations 


A  somewJiat  more  unified  approach  to  generating  highly  stretched  triangulations  involves  the  use 
of  min-max  triangulations,  f.e.,  triangulations  which  minimize  the  maximum  angles  produced  in 
the  elements.  While  Delaunay  triangulations  were  written-ofT  as  unsuitable  for  stretched  mesh 
generation,  the  desire  to  avoid  obtuse  or  large  angle  triangles  in  such  meshes  implies  the  feasibility 
of  using  min-max  triangulations  for  generating  stretched  meshes. 

However,  connectivity  strategies  alone  cannot  guarantee  high  quality  meshes,  and  thus  com¬ 
patible  point-placement  strategies  must  be  developed.  The  advancing-front  min-max  triangulation 
algorithm  of  Marcum  [43]  makes  use  of  two  distinct  point  placement  strategies,  one  for  highly 
stretched  regions  of  the  mesh,  and  another  for  isotropic  regions.  In  the  stretched  mesh  regions, 
the  point  placement  strategy  resembles  that  described  in  the  advancing-layers  method,  ne.,  points 
are  placed  along  smoothed  boundary  normal  segments.  In  the  isotropic  regions,  points  are  placed 
according  to  the  same  rules  described  for  the  advancing-front  Delaunay  triangulation  algorithm 
of  section  3.5,  making  use  of  a  background  element-size  distribution  function.  The  mesh  gener¬ 
ation  process  begins  by  generating  a  coarse  triangulation  which  covers  the  entire  domain.  The 
domain  boundaries  are  then  discretized  according  to  the  background  function,  and  these  boundary 
points  are  inserted  into  the  existing  triaiigulation.  Once  all  boundary  points  have  been  inserted, 
the  boundary  integrity  is  recovered.  The  advancing-front  min-max  procedure  is  then  initiated  by 
creating  new^  points  using  the  ap])ropriate  point-placement  strategy  and  inserting  each  point  into 
the  existing  triangulation  using  the  min-max  variant  of  the  Green-Sibson  algorithm  forced 

triangulation  followed  by  edge-swapping).  The  unifying  feature  of  this  strategy  is  that  the  inser¬ 
tion  and  triaiigulation  of  new  points  is  identical  in  stretched  and  unstretched  regions.  The  point 
placement  methods  are  however  different  in  these  tw^o  regions,  and  the  method  retains  a  somewhat 
hybrid  character.  An  example  of  a  stretched  unstructured  mesh  generated  by  this  technique  is 
shown  in  Figure  34. 


Figure  34:  Stretched  unstructured  mesh  produced  by  advancing-front  min-max  trian¬ 
gulation  procedure  for  familiar  two-dimensional  geometry,  with  illustration  of  stretched 
triangulation  over  the  coast  of  Belgium.  (Reproduced  from  [43]  wdth  permission). 


The  method  has  also  been  demonstrated  in  three  dimensions,  as  shown  in  Figure  35.  In  three 
dimensions,  each  time  a  new'  point  is  inserted,  the  mesh  must  first  be  sw'apped  to  a  local  Delaunay 


40 


configuration,  and  then  to  the  min-max  configuration,  in  order  to  attain  a  rnore  globally  optimum 
final  min-max  configuration,  as  described  in  section  3.5. 

An  alternate  point-placement  strategy  for  min-max  triangulations  which  does  not  rely  on  an 
advancing-front  is  given  by  Barth  [65].  One  of  the  drawbacks  of  the  reliance  on  min-max  tri¬ 
angulations  is  the  possibility  of  the  min-max  criterion  resulting  in  undesirable  connectivities  for 
seemingly  well  distributed  vertices.  For  highly  stretched  meshes,  the  min-max  triangulation  may 
be  very  sensitive  to  the  placement  of  vertices,  with  smaU  displacements  leading  to  substantially 
different  connectivities.  The  precise  placement  of  grid  points  in  the  highly  stretched  regions  is  thus 
much  more  critical  than  in  “forced”  triangulation  methods  such  as  the  advancing-layers  method. 


Figure  35:  Section  of  an  advancing-front  min-max  triangulation  grid  in  three  dimensions 
about  an  F-18  aircraft  configuration.  (Reproduced  from  [43]  with  permission). 
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